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ABSTRACT 


rJ 

ine  goal  of  this  thesis  was  to  develop  and  validate  a  set  of 
microcomputer  programs  using,  in  part,  a  previously  written 
finite-element  algorithm  to  analyze  the  perturbation  of  an  incident 
electromagnetic  field  as  it  penetrates  a  missile  radome.  An 
interactive  program  was  developed  to  design  the  radome  structure 
using  aerodynamic  shaping  functions  and  provide  structure  and 
material  files  as  input  to  the  finite-element  code.  A  second  program 
was  developed  to  use  the  spherical  harmonic  expansion  coefficients 
provided  by  the  finite-element  code  to  assemble  the  electromagnetic 
fields  within  the  radome  core  and  to  display  the  fields  that  appear 
across  a  planar  antenna  with  three-dimensional  graphics  for  any 
orientation  of  the  antenna.  Algorithms  were  included  which  compare 
the  computed  field  components  to  the  theoretical  incident  plane 
wave  for  each  stipulated  angle  of  incidence,  in  order  to  determine 
the  perturbation  due  to  the  presence  of  the  radome.  Validation  of 
the  computational  method  was  attempted  u  nalyzing  the 

perturbation  indicated  for  an  ideal  radome  with  relative  permittivity 
of  unity.  The  attempted  validation  showed  phase  errors  in  the 
computed  fields  which  were  minimal  for  axial  incidence,  but  became 
significant  for  highly  canted  incidence. 
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I.  INTRODUCTION 


Radome  design  is  a  composite  discipline  evolved  from 
aeronautical,  electrical,  mechanical  and  materials  engineering.  The 
word  radome  was  coined  during  World  War  II  as  a  composite  of  two 
words,  radar  and  dome  [Ref.  l:p.  ix].  The  term  radome  originally 
described  the  family  of  dome-shaped  structures  designed  and  built 
to  house  and  protect  the  radar  antennas  on  aircraft,  but  now  the 
term  is  applied  to  any  structure  that  protects  a  radio  frequency 

antenna  system.  The  structure  is  designed  to  provide  maximum 
protection  to  the  antenna  with  minimum  interference  to  the 
antenna's  electrical  performance. 

Radomes  on  missiles  are  presently  designed  to  meet  many 
diverse  specifications.  These  specification  requirements  can  be 

divided  into  two  categories,  electrical  and  aero-mechanical.  Electrical 
design  problems  include  those  that  have  to  do  with  the  radome’s 
ability  to  act  as  a  window,  ideally  transparent,  to  the  electromagnetic 
energy  that  must  pass  through  the  radome  to  be  used  by  the 

missile's  antenna  system.  This  energy  is  of  a  specific  frequency 
range  and  intensity  that  is  within  the  operating  parameters  of  the 

antenna  system.  The  magnitude  and  form  of  the  electromagnetic 
energy  will  be  dependent  on  range  from  the  source,  polarization,  and 
incident  angle.  Aero-mechanical  design  problems  address  not  only 
the  design  of  the  structure  of  the  radome,  but  also  the  ability  of  the 


1 


radome  to  withstand  environmental  conditions  and  battle  damage. 
Environmental  conditions  vary  with  scenario,  and  include:  normal 
pressure  force  and  viscosity,  wind  loads,  inertial  loads,  temperature 
variations,  moisture  resistance,  resistance  to  rain  erosion,  and  other 
diverse  factors.  Design  of  a  radome  for  use  in  battle  may  require 
that  the  aero-mechanical  design  consider  resistance  to  nuclear 
radiation  and  laser  weapons. 

As  in  any  engineering  design,  there  are  tradeoffs  between  the 
various  requirements,  such  that  each  interrelated  component  design 
factor  must  be  weighted  and  considered  in  the  overall  system  design. 
Numerical  models  of  environmental  conditions  assist  in  the  analysis 
of  the  tradeoffs  for  the  aero-mechanical  design.  Numerical  models  to 
aid  in  the  electrical  design  analysis  of  a  radome  to  protect  missile 
antenna  systems  which  operate  at  high  frequencies  are  available, 
such  as  the  ray-tracing  algorithms  described  by  Hirsch  and  Grove  in 
Ref.  2.  When  a  missile  antenna  system  requires  operation  in  the 
resonance  region  of  frequencies  for  the  radome  parameters,  a 
numerical  model  is  required  that  can  solve  the  radome  interaction 
problem. 

The  subject  of  this  thesis  is  the  development  of  a  generalized, 
interactive  set  of  programs  to  be  used  in  the  analysis  of  the 
interaction  between  an  actual  or  proposed  radome  structure  and 
incident  electromagnetic  fields  at  the  lower  frequencies  of  interest, 
where  strong  resonance  characteristics  exist.  In  order  to  analyze  a 
specified  radome,  a  numerical  model  is  constructed  using  the 
interactive  program,  RADSHP.  This  radome  generation  program 
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provides  geometric  shaping  functions  to  allow  the  user  to  construct  a 
multilayer  radome  model  of  specified  materials.  Chapter  II  provides 
a  discussion  of  radome  design  methodology,  the  aerodynamic  shaping 
functions  available  in  this  package,  radome  wall  constructions,  and 
material  specification. 

This  thesis  work  is  an  application  of  the  programs  developed 
for  use  on  a  personal  computer  by  Morgan  [Ref.  3]  and  Connolly 
[Ref.  4]  to  calculate  electromagnetic  scattering  by  inhomogeneous 
dielectric  and/or  magnetic  bodies  of  revolution.  As  discussed  in 

Chapter  II,  radomes  are  manufactured  using  dielectric  materials  in 
order  to  provide  a  radio  frequency  window  to  the  antenna  housed  in 
the  radome’s  interior.  An  axisymmetric  missile  radome  structure 

meets  the  physical  structure  requirement  for  rotational  symmetry. 
The  numerical  analysis  of  the  interaction  of  the  specified  radome 
structure  and  the  incident  electromagnetic  energy  is  performed  by 

the  code  developed  by  Morgan,  EMCAD.  Modifications  to  this  code 
resulted  in  the  program  EMRAD,  which  provides  data  file  generation 
to  store  the  modal  expansion  coefficients  which  are  used  by  the 

program  CORFLD  to  determine  the  electromagnetic  fields  within  the 
radome  interior.  In  Chapter  III,  the  theory  required  for  the 
formulation  of  the  computed  interior  fields  of  the  radome  by  the 
program  CORFLD  is  discussed. 

The  programs  developed  by  Morgan  use  an  optimized 
variational  finite-element  algorithm  in  conjunction  with  a 
tri-regional  unimoment  method  to  provide  scattering  solutions  for 
each  of  multiple  incident  fields  impinging  on  the  specified  structure. 
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The  programs  represent  a  significant  improvement  in  computational 
speed  and  versatility  over  computer  codes  which  find  the  solution  to 
the  electromagnetic  resonance  behavior  by  use  of  surface  current 
integral  equation  formulations.  An  overview  of  the  finite  element 
scattering  algorithm,  along  with  the  concepts  embodied  in  the 
coupled  azimuthal  potentials  and  the  unimoment  method  are 
included  in  Appendix  A.  Further  information  is  available  in 
references  3,  4,  5,  6,  and  7. 

Chapter  IV  discusses  the  flow  of  data  between  the  programs 
which  have  been  developed.  The  general  use  of  the  RADSHP 
program  to  construct  the  radome  data  files  is  discussed,  as  well  as 
the  data  files  generated  and  their  uses.  The  flow  of  information  into 
and  out  of  the  radome  analysis  program  EMRAD  is  analyzed.  The 
execution  and  operation  of  CORFLD  is  discussed. 

In  Chapter  V,  the  validation  of  the  interior  field  computation  is 
discussed,  wherein  the  numerically  generated  fields  are  compared  to 
the  incident  fields  for  a  radome  constructed  of  air.  The  results 
produced  for  radomes  constructed  of  materials  presently  used  in  the 
radome  industry  are  presented  and  discussed. 

The  conclusions  in  Chapter  VI  will  summarize  the  work 
previously  presented,  as  well  as  provide  recommendations  for  future 
modifications  to  the  interactive  package  of  programs  being  used. 

Additional  Appendices  are  included  to  provide  copies  of  the 
source  code  for  the  computer  programs,  RADSHP  and  CORFLD,  which 
were  written  as  part  of  this  thesis  effort.  Source  listings  for 
additional  programs,  such  as  EMCADIN,  EMESH,  and  EMRAD  are 
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available  upon  request  from  Prof.  Michael  A.  Morgan.  His  address  is 
given  in  Appendix  D,  Software  Sources. 


5 


II. 


A.  RADOME  DESIGN  METHODOLOGY 

The  configuration  of  a  radome  for  high-speed  missiles  must 
combine  high  ae:  adynamic  performance  with  a  highly  accurate 
transmission  of  electromagnetic  energy  through  the  radome  to  the 
missile  seeker.  As  described  in  Ref.  2,  the  methodology  for 
establishing  radome  performance  requirements  is  an  iterative 
process.  Figure  2.1  shows  the  process  as  a  set  of  design  procedure 
loops  that  begin  with  a  statement  of  desired  performance,  and 
include  analytical  studies,  hardware  development,  and  hardware 
evaluation.  These  result  in  a  set  of  optimized  radome  requirements, 
which  is  the  blueprint  from  which  a  radome  can  be  constructed. 
Computer  software  tools  help  to  maximize  the  benefit  of  the  work 
done  in  analytical  studies,  thereby  minimizing  the  number  of 
iterations  required  for  the  hardware  development  and  evaluation. 
The  programs  described  in  this  thesis  are  designed  to  analyze  the 
interaction  between  the  aerodynamic  considerations  in  radome 
design  and  the  electrical  performance  of  the  radome. 

The  aerodynamic  considerations  in  the  design  of  a  radome 
structure  are  three-fold:  shape-generated  drag,  thermal  stress  from 
aerodynamic  heating,  and  structural  performance.  These  three 
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parameters  are  not  independent  of  each  other,  and  the  electrical 
performance  of  the  radome  is  a  strong  function  of  these  aerodynamic 
factors.  The  outer  shape  of  the  radome  has  a  large  impact  on  how 
the  missile  will  fly  through  the  air.  It  is  this  outer  conformation  that 
will,  along  with  speed  of  travel  and  angle  of  attack,  determine  the 
forces  that  act  on  the  radome  structure  during  flight.  The  forces 
generated  during  flight  determine  the  design  requirements  on  the 
radome  structure  for  material  strength  and  thermal  stress. 
Fortunately,  the  shapes  that  result  in  the  best  drag  performance  are 
those  that  can  also  be  used  to  develop  radomes  that  match  the 
structural  requirements.  Unfortunately,  the  shapes  that  best  meet 
the  drag  requirements  are  not  necessarily  the  best  for  maximizing 
the  electrical  performance.  The  programs  developed  in  this  thesis 
include  the  computer  aided  design  (CAD)  functions  to  analyze  several 
shaping  functions  that  aerodynamicists  have  developed  to  minimize 
the  aerodynamic  drag  forces  with  respect  to  radome  structure 
parameters. 

B.  RADOME  SHAPING  FUNCTIONS 

An  ogive  shaping  function  is  one  of  the  most  common  used  in 
radome  designs.  The  ogive  shape  is  popular  because  it  is  easily 
fabricated  and  has  relatively  acceptable  aerodynamic  performance  at 
high  speeds.  It  is  defined  {Ref.  l:p.  48]  as  a  segment  of  the  arc  of  a 
circle,  whose  radius,  F},  is  greater  than  the  radius  of  the  radome 
base.  A  pointed  ogive  is  generated  when  the  circle's  radius,  Fj,  is 
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large  enough  so  that  the  arc  reaches  the  radome’s  axis  of  symmetry. 
The  ogive  shaping  function  is  described  by 


x=  y  Fj2  -  ((H-z)-F2  )2  -  F3  ,  (2.1) 

where  F2  and  F3  determine  the  position  of  the  circular  shaping  arc,  as 
illustrated  in  Fig.  2.2.  The  ogive  shaping  function  degenerates  to  a 
hemisphere  when  the  shaping  factors,  F2  and  F3,  are  defined  as 
F2  =  Fj  and  F3  =  0. 

A  tangent  ogive  shaping  function  generates  a  radome  shape 
such  that  the  surface  of  the  radome  is  tangent  at  the  base  of  the 
radome,  providing  a  smooth  transition  to  the  missile  body  [Ref.  8].  A 
tangent  ogive  shaping  function  is  determined  such  that 


where 


(2.4) 


A  shaping  function  that  provides  maximum  volume-to-drag 


ratio  at  high  flight  speeds  is  the  von  Karman  function  [Ref.  1  :p.  48]  , 


j  sin(2Fz ) 


1 

2 


(2.5) 
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Radome  shape  -  nomenclature 


where 


Fz  =  cos  1  1 


9  (H~Z> 

H 


(2.6) 


Radomes  described  by  this  von  Karman  shaping  function  have  a 
limited  application  due  to  fabrication  difficulties  related  to  the 
complex  shape. 

In  order  to  minimize  drag  for  a  given  height-to-radius  ratio,  a 
Newtonian  contour  is  used.  This  contour  is  approximated  by  the 
power  series  radome  shaping  function  [Ref.  l:p.  49], 


x  =  R 


(H~z) 

H 


^6 


(2.7) 


The  exponent,  F6,  provides  for  shaping.  When  F6  =  1.0,  the  radome 
shape  generated  is  conical.  Aerodynamicists  normally  use  F6  =  0.75 

to  approximate  the  Newtonian  contour. 

A  shaping  function  analyzed  by  Hirsch  and  Grove  [Ref.  2]  uses  a 
parabolic  basis  function  for  the  radome  shape. 


x2  =  F7-z  .  (2.8) 

The  rounded  effect  developed  at  the  tip  of  the  radome  by  the 
parabolic  shape  function  can  be  eliminated  by  the  creation  of  a 
mathematical  "deadzone"  for  points  within  a  region  about  the 
radome's  apex.  The  parabolic  shaping  coefficient,  F7  ,  determines  the 
extent  of  the  "dead  zone”  and  thereby  controls  the  slope  of  the 
parabolic  surface.  The  radome  shape  is  determined  in  terms  of  the 

shaping  coefficient,  then  scaled  to  the  radome  size  parameters.  The 
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unsealed  shaping  function. 


x  =  Jf7  -z  ,  (2.9) 

is  used  to  generate  a  parabolic  radome  shape  with  a  height  of  F7  .  In 
order  to  eliminate  the  parabolic  rounding  of  the  apex,  the  "deadzone" 
is  established  for  any  points  for  the  unsealed  z  value  within  1.0  of  F7. 
Within  this  deadzone,  z  >  (  F7  -  1  ),  the  radome  shape  is  determined 
at  a  line  of  constant  slope,  converging  at  a  height  of  z  =  F7  +  1.0.  This 
results  in  the  marrying  of  a  pointed  apex  onto  the  parabolic  radome 
shape,  with  the  magnitude  of  the  shaping  coefficient,  F7  ,  determining 

the  relative  magnitude  of  the  apex  to  the  total  radome  height. 

This  project  contains  coding  in  RADSHP  to  generate  radomes  of 
multiple  layers  for  the  shaping  functions  described:  the  general 
ogive,  the  tangent  ogive,  the  von  Karman  shape,  the  power  series 
radome,  and  the  parabolic  radome  with  pointed  apex.  RADSHP  can 
generate  a  complex  radome  construction  model  using  different 
shaping  functions  for  each  layer.  RADSHP  generates  a  material  data 
file  appropriate  for  input  to  the  numerical  analysis  program,  EMRAD, 
and  a  structure  data  file  that  can  be  viewed  and  modified  using 
CURVE  DIGITIZER  (Appendix  2)  ,  then  translated  by  EMCADIN  into  an 
input  structure  data  file  for  EMRAD. 

C.  RADOME  WALL  CONSTRUCTION 

If  a  material  existed  which  was  mechanically  strong  enough  to 

sustain  the  aerodynamic  forces  that  the  radome  experiences  during 
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flight,  and  that  same  material  had  the  electromagnetic  properties  of  a 
vacuum,  then  the  radome  could  be  built  to  ideally  meet  all 
requirements.  Such  a  material  does  not  exist,  and,  in  general,  as  the 
requirement  for  material  strength  increases,  the  electromagnetic 
properties  tend  to  strongly  deviate  from  the  properties  of  the  ideal 
vacuum.  The  transmissivity  of  the  radome  tends  to  decrease  as  the 
dielectric  constants  of  the  construction  materials  increase.  In  order 
to  meet  strength  requirements,  while  minimizing  the  perturbation 
caused  to  the  electromagnetic  field,  layer  construction  of  radomes 
has  been  developed. 

Lamination  construction  methods  allow  for  an  alternating 
combination  of  thin  layers  of  high  density,  high  strength  materials  of 
high  dielectric  constants,  with  thicker  layers  of  low  density,  lower 
strength  materials  of  low  dielectric  constants.  Electromagnetic  design 
of  laminated  radomes  is  more  complex  than  the  design  of  monolithic 
radomes  in  that  there  is  a  greater  number  of  boundary  conditions  to 
be  considered  in  establishing  the  lamination  layer  thicknesses.  The 
EMRAD  program  developed  by  Morgan  [Ref.  3]  currently  handles 
layered  dielectrics  with  up  to  five  layers. 

Figure  2.3,  adapted  from  the  Radome  Engineering  Handbook 
[Ref.  l:p.  46],  shows  several  types  of  laminated  wall  construction  that 
exhibit  the  alternation  of  layers  of  low  dielectric  constant  and  high 
dielectric  constant  materials.  The  two-ply  sandwich  is  the  simplest 
form  of  lamination  construction  with  a  thin  skin  of  strong,  dense 
material  of  high  dielectric  constant,  supported  internally  by  a  thicker 
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layer  of  porous  material  of  low  dielectric  constant.  This  design  is 
especially  attractive  at  low  frequencies  as  the  lateral  dimension  of 
the  thin  layer  of  dense  material  is  only  a  fraction  of  the  longer 

wavelength. 

The  A-sandwich  and  B-sandwich  are  two  variations  of  a 
stronger  three-layer  construction.  Although  the  A-sandwich  is 
superior  electrically  to  the  B-sandwich,  while  the  placement  of  the 
high  density  material  in  the  center  of  the  layer  construction  will 

strengthen  the  radome  wall  by  making  it  stiffer,  the  porous,  low 
density  material  on  the  exterior  of  the  radome  wall  will  not  sustain 
the  aerodynamic  forces  during  flight. 

The  B-sandwich's  inner  and  outer  skin  of  high  density  material 
and  the  core  of  low  dielectric  material  provide  stiffness  and  strength 
to  the  radome  structure  while  providing  a  dense  outer  skin.  This 
construction  has  a  higher  strength-to-weight  ratio  than  a  monolithic 
construction.  The  core  thickness  is  determined  by  the  strength 
requirements  and  provides  spacing  between  the  inner  and  outer 

skin.  Optimal  electrical  spacing  between  the  skins  would  by 

frequency  dependent  at  a  core  thickness  of  a  quarter  wavelength 
[Ref  l:p.  45]. 

The  C-sandwich  and  multilayer  radome  wall  constructions 
provide  for  very  high  strength  levels  which  may  be  required  at  high 
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Figure  2.3  Laminated  radome  wall  construction 


operating  speeds  or  for  use  in  inclement  weather.  The  layer  material 
may  consist  of  plastic  laminates  and  ceramics.  Multilayered 
construction  is  attractive  where  extreme  structural  rigidity  or 
broadband  capabilities  are  primary  requirements.  [Ref.  l:p.  46] 

The  radome  interacts  with  the  electromagnetic  field  as  a  lens, 
modifying  the  field  as  it  passes  through  the  radome  structure.  At 
higher  frequencies,  this  interaction  can  predicted  by  the  use  of 
ray-tracing  techniques,  similar  to  optical  lens  analysis  [Refs.  2  and 
8].  Studies  have  been  conducted  at  higher  frequencies  to  design  a 
layer  of  dielectric  material  within  the  radome’s  interior  structure  to 
act  as  a  correcting  lens  [Ref.  8]  as  shown  in  Fig.  2.4  .  This  package  of 
programs  will  allow  similar  analysis  of  lens  construction  at  lower 
frequencies,  where  the  radome  and  lens  parameters  are  within  the 
resonance  region  of  the  electromagnetic  fields. 

D.  MATERIALS  FOR  RADOME  CONSTRUCTION 

The  choice  of  material  for  construction  of  a  radome  is  another 
engineering  design  decision.  As  material  science  engineers  continue 
to  develop  new  materials,  particularly  low-weight,  high-strength 
ceramics,  the  range  of  materials  available  to  the  radome  designer 
increases.  The  general  properties  of  the  three  groups  of  materials 
which  are  in  common  use  are  listed  in  Ref.  9,  and  is  used  here  as  a 
source  for  material  selection.  The  groups  of  materials  discussed  are 
fiber-reinforced  organic  resins,  refractory  oxides,  and  glass  ceramics. 
Physical  properties  are  listed  in  Table  2.1. 


Relative 

permittivity 

K ) 


Loss  tangent 
( tan  5 ) 


Density 


Ultimate 

Tensile 

Strength 

(UTS) 

Compressive 

Strength 


TABLE  2.1 

RADOME  MATERIAL  PROPERTIES 
FIBER 

REINFORCED  GLASS  REFRACTORY 

RESINS  CERAMICS  OXIDES 


4-6  5.5  8-9 

0.01  0.003  0.002 


N/A 


160  lb/ft3  220  lb/ft3 


40  kpsi  7  kpsi  27  kpsi 


20-40  kpsi  35  kpsi  280  kpsi 


Maximum 

Temperature 


260  °C 


500  °C 


1700  °C 


In  fiber-reinforced  organic  resins,  a  variety  of  resin  types  are 
used  to  impregnate  fiberglass.  Organic  resins,  such  as  phenolic, 
polyester,  and  epoxy  are  used  with  the  fiber  to  construct  the  radome. 
They  are  then  thermally  treated  to  set  the  structure.  Several 
manufacturing  techniques  are  discussed  in  Ref.  1.  Control  of  layer 

thickness  is  difficult  to  maintain  and  requires  finishing  to  avoid 

aberration.  Radomes  constructed  of  these  materials  have  poor 

resistance  to  rain  erosion,  but  are  relatively  cheap  and  easy  to 

manufacture,  even  in  large  sizes.  Maximum  temperatures  that  the 
material  can  sustain  limit  use  of  these  materials  to  missiles  with 
speed  of  less  than  Mach  1.5  to  Mach  2.0. 

Typical  glass  ceramics  used  in  radome  construction  are 
Pyrocream  and  glass-bonded  mica.  Compared  to  organic  radome 
materials,  ceramic  materials  have  much  higher  resistance  to  rain 
erosion  at  supersonic  speeds.  Ceramics  are  brittle  materials  and  the 
controls  needed  during  manufacturing  to  ensure  structural  integrity 
of  the  radomes  makes  them  expensive  to  build  [Ref.  l:p.  242], 

The  refractory  oxide  used  in  radome  construction  is  usually 
alumina  (A1203),  but  magnesia  (MgO)  and  fused  silica  (Si02)  are  also 
used.  Refractory  oxide  manufacturing  is  also  expensive  due  to  the 
high  temperatures  used  in  forming  these  materials  to  the  desired 
radome  shape.  These  materials  have  high  dielectric  constants,  that 
can  undergo  a  change  of  up  to  ten  percent  between  room 
temperature  and  1700  °C.  Some  of  the  methods  by  which  refractory 
oxide  radomes  are  manufactured  are  slip  casting,  isostatic  pressing. 


and  flame  spraying  [Ref.  l:pp.  298-313].  Oxides  have  a  high 
resistance  to  rain  erosion  and  wall  thickness  can  be  controlled  to 
within  0.002  inches,  allowing  for  extreme  control  of  aberration. 
Although  the  dielectric  constant  of  these  materials  is  high,  they 
provide  the  highest  degree  of  thermal  protection  and  place  the  least 
restriction  on  missile  speed. 
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III.  ELECTROMAGNETIC  FIELD  THEORY 


A.  BACKGROUND 

The  analysis  of  the  interaction  between  a  missile's  radome  and 
the  electromagnetic  signal  that  is  radiated  or  received  by  that 
missile's  radar  is  a  problem  similar  to  those  undergoing  intensive 
investigation  in  electromagnetic  interaction  with  other  physical 
medium.  Ongoing  study  of  medium  interactions  include  other 
military  weapons  and  weapons  platforms,  forested  areas,  bodies  of 
water,  rain,  and  biological  tissues  [Refs.  3,  5,  and  10],  The  study  of 
electromagnetic  wave  penetration  through  missile  radomes  allows 
for  a  special  approach  to  the  solution  of  the  problem  for  the  common 
case  of  rotational  symmetry  in  the  construction  of  the  radome.  In 
the  high-frequency  regime,  the  solution  may  be  found  as  described 
by  Hirsch  and  Grove  [Ref.  2],  with  application  of  an  asymptotic 
approach  using  the  geometrical  theory  of  diffraction  (GTD),  which 
incorporates  ray-tracing.  When  lower  frequencies  are  considered, 
and  the  interaction  problem  is  in  the  resonance  region,  a  complete 
solution  of  Maxwell’s  equations  is  required.  The  asymptotic  GTD 
method  yields  unsatisfactory  results. 

The  solution  of  Maxwell's  equations  for  the  radome  used  within 
this  project  is  divided  into  three  regions  as  illustrated  in  Fig.  3.1.  The 
three  regions  are  the  homogenous  spherical  inner  radome  core 
(designated  Region  I);  the  transition  region,  where  a  semi-annular 
conformal  mesh  is  generated  for  field  analysis  (designated  Region  II); 
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and  an  unbounded  exterior  scattering  region  (designated  Region  III). 
The  total  solutions  for  the  fields  in  the  core.  Region  I,  and  the 
scattered  fields  in  the  exterior  region,  Region  III,  are  produced  via 
separate  spherical  harmonic  expansions.  Using  the  unimoment 
method  (Appendix  A),  the  modes  in  these  expansions  are  applied  as 
boundary  conditions  on  the  inner  and  outer  spherical  surfaces 
bounding  the  mesh  and  the  resultant  fields  for  each  such  mode  are 
found  using  the  finite  element  method.  A  similar  analysis  is 

performed  for  the  incident  fields  on  the  outer  surface. 

Using  superposition,  the  finite  element  solutions  are  then 
combined  using  the  same  unknown  coefficients  as  in  the  original 
expressions.  The  combined  finite  element  solution  is  equated,  in  the 
least  squares  sense,  to  the  original  expressions  along  the  inner  and 
outer  surfaces.  Enforcement  of  this  equality  allows  solution  of  the 
expansion  coefficients  for  the  core  and  exterior  regions. 

The  exterior  expansion  coefficients  are  used  by  EMRAD  to 
assemble  the  scattered  fields  and  bistatic  radar  cross  sections  (see 
Refs.  3  and  4).  The  core  region  spherical  harmonic  expansion 

coefficients  are  used  to  construct  the  fields  that  penetrate  the 
radome  to  illuminate  the  enclosed  antenna.  The  details  of  how  these 
interior  region  fields  are  generated,  using  the  core  region  spherical 
harmonics,  is  the  topic  of  the  next  section.  An  overview  of  the 

theory  used  by  Morgan  in  EMRAD  is  included  in  Appendix  A. 

Further  information  is  available  in  Refs.  3,  4,  5,  6,  and  7. 
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B.  GENERATION  OF  THE  CORE  REGION  FIELDS 

In  this  section  we  will  consider  the  construction  of  the  core 
region  fields  using  the  expansion  coefficients  provided  by  EMRAD. 
These  fields  are  assembled  and  plotted  in  various  ways  over 
user-selected  planar  surfaces  by  CORFLD. 

Within  the  spherical  core,  the  source-free  phasor  Maxwell’s 
curl  equations,  with  suppressed  e  ^  time  variation,  are 

-VxE=jco|iH  (3.1) 

and 

VxH=jtoeE  ,  (3.2) 

where  p  and  e  are  the  respective  permeability  and  permittivity.  The 
p  and  e  will  usually  be  those  of  free-space  for  the  case  of  the  core 
region  in  a  radome. 

In  a  homogenous,  source-free  spherical  region,  such  as  the 
radome  core,  it  can  be  shown  [Ref.  11]  that  the  vector 
electromagnetic  field  can  be  represented  by  the  combination  of  two 
types  of  spherical  harmonic  fields,  transverse  magnetic  (TM)  and 
transverse  electric  (TE).  Using  the  spherical  coordinates  of  Fig.  3.2, 
these  TM  and  TE  fields  have  the  respective  properties  that  Hr  =  0  and 
Er  =  0.  Furthermore,  the  TM  and  TE  vector  fields  can  be  generated 
from  radially  directed  vector  potentials,  such  that  for  the  TM  fields, 

A  =  Ar  r  (3.3) 

and 
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F  =  Fr  r  , 

where  the  radial  components  satisfy 


(3.4) 


'  Ar 


+k2 


V 


0 


(3.5) 


and 


vr ; 


+k2 


vr  ; 


=  o  . 


(3.6) 


These  equations  are  in  the  form  of  the  complex  scalar 
Helmholtz  equation, 

V2vP  +k2vF  =  0  .  (3.7 ) 

In  spherical  coordinates,  the  Helmholtz  equation  becomes  [Ref.  11], 


i-i. 

t2  0r 


an* 

3r  J 


1  ^  1  — +k2  ^=0 

l+r2  sine  ^01  ae  /  r2sin20  34>2 


(3.8) 


To  find  a  solution  to  this  differential  equation,  the  method  of 
separation  of  variables  is  used.  Letting 

S’  =  R(r)  8(8)  <K»)  ,  (3'9) 

the  separation  of  variables  procedure  yields  the  system  of  ordinary 

differential  equations: 

dJr2^j,{k2r2-n(n+l)|R  =  0  ,  (3.10) 
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sin20 


6  =  0  , 


(3.1  1) 


and 


d<D2  -> 

— — +m2 O  =  0  .  (3.12) 

d<|>2 

The  R  equation  is  related  to  Bessel's  equation,  and  its  solutions  are  in 
the  form  of  spherical  Bessel  functions.  Spherical  Bessel  functions,  bn  , 
are  related  to  ordinary  Bessel  functions,  Bn,  as  indicated. 


R(r)  =  bn(kr) 


(3.13) 


The  0  equation  is  related  to  Legendre’s  equation  and  its  solutions  are 
the  associated  Legendre  functions, 


0(0)  =  L™(cos0)  ->  F™(cos0)  ,  Q£(cos0)  ,  (3.14) 

where  Pnm(cos  0)  are  the  associated  Legendre  functions  of  the  first 
kind  and  Qnm(cos  0)  are  the  associated  Legendre  functions  of  the 

second  kind.  The  associated  Legendre  functions  of  the  second  kind 
tend  toward  infinity  at  points  along  the  z-axis  and  must  be  excluded 
for  solutions  which  include  this  region. 
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The  O  equation  has  harmonic  sinusoidal  type  solutions 
$(<}>)  =  h(m<>)  — >  cos(m<{>),  sin(m<}>),  e±jm<s 
=  q  cos(m<}))+ C2sin(m<{>) 

=  C3e+jm<>+ C4e~■>m<^,  .  (3.15) 

Together,  these  separated  solutions  form  the  product  solutions 

to  the  Helmholtz  equation  as 

¥ro,n  =  bn(kr)  Pnm(cos6)  h(m 0)  .  (3.16) 

Equations  of  this  type  form  the  elementary  spherical  harmonic  wave 
functions.  The  general  solution  to  the  Helmholtz  equation  is  formed 
by  summation  of  linear  combinations  of  these  elementary  functions, 
one  for  the  TM  fields  and  one  for  the  TE  fields.  Each  of  these 

solutions  will  be  of  the  form 

nr,e,0)  =  XX  MM  P^(cos0)  h (m <J>)  .  (3.17) 

m  n 

The  summation  over  all  possible  m  and  n  is  limited  by  the 

nature  of  the  Legendre  polynomial  of  the  first  kind  in  that  Pnm  =  0 

for  n  <  m.  The  summation  is  also  limited  by  the  restriction  that  n 
will  not  equal  zero  when  m  is  zero.  For  both  n  and  m  equal  to 
zero,  then  the  field  evaluated  is  zero,  therefore  this  solution  is  also 
excluded.  Then  m  is  zero  and  the  summation  over  n  will  start 
with  n  equal  to  one.  This  leads  to  the  scalar  field  expansion 
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oo  oo 

'P(r,0,<J>)=£  X(cmnCosrn<{)  +dmn  sinm$)  P™(cos0)  bn(kr).  (3.18) 

m=0  n=m 
n*0 

As  the  Helmholtz  equation  under  consideration  is  given  as 

V2*F  +k2vP  =  0,  (3.1  9) 

then  for  the  TM  case 

V  =  ,  (3.20) 

r 

and  for  the  TE  case 

'F  =  -5-  .  (3.2  1) 

r 

The  spherical  harmonic  solutions  of  the  Helmholtz  equations  for 
these  sets  of  modes  are  then  of  the  form 

Ar 

■  =  =  l(kr)  bn (k  r) ]  Pnm(cos0)  h(m0)  .  (3.22) 

Fr  . 

In  this  renormalized  form  of  the  solution  to  the  scalar  Helmholtz 
equation,  the  Riccati  (or  Schelkunoff)  spherical  Bessel  functions  are 
present.  The  Riccati  spherical  Bessel  functions  relate  to  the  spherical 
Bessel  function  and  normal  Bessel  functions  as  follows: 


Jn(kr)  =  [  (kr)bn  (kr)  ]  = 


B,  ,  (kr)  . 

K) 


(3.2  3) 


Therefore,  in  the  spherical  core  region,  the  solutions  for  Ar  and  Fr  are 
of  the  form 
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=  X  £(cmnc°sm<)>  +dmn  sinm<)>)  P™(cos0)  Jn(kr)  .  (3.24) 

m=0  n=m 
n*0 

When  expressed  in  terms  of  exponential  Fourier  series  in  0,  the 
solutions  are  of  the  form 

=  £  £amn  Pnm(cos0)  J„(kr)  e'm,>  .  (3.25) 

m=-~  n=m 
n*0 

The  coefficients  within  the  two  summations  are  determined  by  the 
solution  of  the  scattering  problem  provided  by  EMRAD. 

The  spherical  coordinate  electric  and  magnetic  field 
components  can  now  be  found  for  the  TE  and  TM  cases.  For  the  TM 
case,  where 

H™  =  Vx(A,  r)  (3.26) 

and 

E™  =  J_VxVx(Ar  r)  .  (3.27) 

J(OE 

the  field  components  are 


E  -  -L  a  +k2  A 

'  jtOE  [  ar2  ^  +k  Ar  J 

(3.28) 

_  1  1  a2Ar 

™  joe  r  3r80 

(3.29) 

c  1  1  32A, 

^  jo)E  rsin0  8rd0 

(3.3  0) 
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(3.3  1) 


(3.3  2) 
(3.3  3) 

(3.34) 

(3.3  5) 

(3.3  6) 
(3.37) 

(3.3  8) 
(3.3  9) 

(3.40) 

(3.4  1) 


The  total  fields  are  the 
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The  solution  for  the  coefficients  in  the  Ar  and  Fr  expansions  is 

performed  by  the  finite  element  and  unimoment  solution  method 
explained  in  Appendix  A.  The  finite  element  method  is  used  by 
EMRAD  and  applied  to  Region  II  between  the  inner  core  and  outer 
exterior  region.  In  this  finite  element  region  there  is  the 
inhomogeneous  axially  symmetric  structure  of  the  radome. 

EMRAD  generates  solutions  for  incident  plane  waves  with 
Poynting  vectors  at  incident  angles  specified  by  the  program  user, 
such  that  the  Poynting  vector's  direction  is  specified  by  the  unit 
propagation  vector  p  =  sin  a  x  +  cos  a  z.  As  any  plane  wave  can  be 
decomposed  into  the  sum  of  two  orthogonal,  linearly  polarized  plane 
waves,  EMRAD  evaluates  two  plane  waves  for  each  incident  angle,  a 
TxM  plane  wave  and  a  TE  plane  wave,  as  shown  in  Fig.  3.3.  The  TM 
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and  TE  plane  waves  are  orthogonal  to  the  Poynting  vector  and  to 
each  other.  The  TE  plane  wave  is  always  in  the  y  direction.  The  TM 
wave  is  in  the  x-z  plane  and  directed  as  determined  by  the 
expression  cos  ax  -  sin  a  z.  For  incident  plane  waves  specified  in 

this  manner,  the  solution  for  the  component  values  of  the  E  fields  in 

the  core  of  the  radome  can  be  determined  from  taylored  equations 

for  the  TM  and  the  TE  incident  plane  waves.  For  the  TM  incident 

wave,  the  core  E -field  component  values  are  determined  as: 
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For  the  TE  incident  wave,  the  core  E  -field  component  values  are 
determined  as: 

Pm(COS0) 


F  sinO  ^  N  Pm(cos 

E{(tM)=  p  ■  "  2  z  Z  bnn(n  +l)Jn(k  r)~— -q- 

V^r  m=0  n=mL  '  " 


n*0 


N 


Eq (r  ,a<t>) Z  *rn  s * °  (m  <t>)  Z 

K0r  m=0 


n=m 

n*0 


Pm(cos0) 


dP 


m 


anJ"(kl)^0-+b"^(kr)W 


(3.5  1) 


(3.52) 


and 
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P  JT"  M  N 

E0(r,e,<j>)=  -J  —  X  ^mCos(m<}>)  £ 

Kor  m=0  n=m 

n*0 

With  the  an  and  bn  coefficients  supplied  by  EMRAD  for  each  pair  of 

plane  waves,  TM  and  TE,  at  each  user-specified  angle,  CORFLD 
generates  the  interior  phasor  field  components  using  the  equations 
above.  The  total  magnitude  of  the  complex  vector  field  at  any  point 
within  the  radome's  core  is  expressed  as  the  square  root  of  the  sum 
of  the  squares  of  the  complex  absolute  values  of  the  field 
components.  This  total  magnitude  is  computed  and  displayed  for 
each  incident  wave  by  CORFLD,  using  the  an  and  bn  coefficients 
supplied  by  EMRAD. 

In  order  to  determine  the  perturbation  of  the  electromagnetic 
field  due  to  the  presence  of  the  radome,  CORFLD  generates  the  values 
of  the  incident  plane  wave,  then  compares  this  to  the  core  region 
field  components  calculated  by  EMRAD  and  CORFLD.  The  formulas 
for  the  undisturbed  incident  TM  and  TE  planes  waves  are 

E™=E0  e-jkoP*r  (xcosa-z sina)  (3.54) 

and 

E11^  e“jkoP’r  y  .  (3.5  5) 

The  magnitude  factor,  E0  =  1.0,  is  set  by  EMRAD  for  both  TM  and  TE 

waves.  The  unit  propagation  vector,  p,  includes  the  direction  of 
propagation  foi  the  incident  waves  and  is  designated  as 

p  =  (x  sina+z  cosa)  .  (3.5  6) 


dP™  -  Pnm(cos0) 
an  Jn(k  r)— ^-+bn  Tn  (k  r)-2— - - 


sin0 


.(3.53) 


35 


This  vector  is  dotted  with  the  position  vector,  r,  for  each  point  on  the 
plane  of  interest,  as  defined  as 

r  =  (xx  +  yy+  zz)  .  (3.57) 

With  the  dot  product  evaluated  in  the  exponential,  the  TM  and  TE 
waves  are  specified  as 

E™=E0  e'jko(xsina+zcosa)  (xcosa-zsina'  (3.5  8) 

and 


Ete_Eq  e~jko(x  s>na  +zcosa)  y  .  (3.5  9) 

The  relative  magnitude  and  phase  of  the  computed  field  in  the 
direction  of  the  theoretical  field  are  calculated  and  displayed.  This 
"scaled  dot  product"  is  defined  as 

Edot  =  .  (3-6  0) 

where  E  is  the  computed  field,  e;  is  a  unit  vector  pointing  in  the 

same  direction  as  either  the  TM  or  TE  incident  field  specified  in 
Equations  (3.58)  and  (3.59),  and 

Ej  =E0  e-jko(xsina+zcosa)  .  (3.61) 

The  remaining  field  term  that  is  computed  and  displayed  by  CORFLD 
is  the  error  term.  The  error  term  is  that  portion  of  the  computed 
field  which  is  perpendicular  to  the  incident  field  vector.  In 

normalized  form  this  becomes. 


E 


error 


E-(E«ei)ei 

E, 


(3.6  2) 
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The  relative  magnitude  of  the  error  term  is  calculated  and  displayed 
by  CORFLD  as  a  measure  of  electromagnetic  disturbance  caused  by 
the  presence  of  the  radome.  This  built-in  theoretical  field 
comparison  is  used  in  Chapter  V  as  a  measure  of  modelling  accuracy 
by  comparing  the  generated  theoretical  fields  to  the  fields  computed 
and  assembled  by  EMRAD  and  CORFLD  for  an  ideal  radome 
constructed  of  air  (er  =  1). 


C.  CANTING  THE  PLANAR  ARRAY 

The  planar  antenna  within  a  radome  can  be  canted  in  any 

direction  allowed  by  the  gimbal  system  that  is  used  with  the 

antenna.  In  order  to  evaluate  the  field  components  across  a  canted 
planar  antenna  within  the  radome  core,  the  coordinates  of  each 
point  on  the  antenna  must  be  converted  from  the  local  planar 
coordinate  system  to  the  global  radome  core  coordinate  system. 

CORFLD  uses  a  transformation  from  local  planar  to  global  core 
coordinates,  assuming  a  two-axis  gimbal  system  where  the  antenna 

is  rotated  about  the  original  z-axis,  as  shown  in  Fig.  3.4,  followed  by  a 
rotation  about  the  new  y-axis,  as  shown  in  Fig.  3.5. 

The  transformation  matrix  is  formed  by  multiplying  two 
matrices,  one  for  each  axis  rotation  [Ref.  12].  The  transformation  for 
the  rotation  about  the  z-axis  is  provided  by  the  matrix 


x 


y 


z 


-k:osP  -sinp  00 

sinp  cosp  00 

00  00  1.0 


(3.6  3) 
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radome 

z-axis 


Figure  3.5  Rotation  about  the  new  y-axis 
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The  transformation  for  the  rotation  about  the  y-axis  is  provided  by 
the  matrix 


V 

cosy 

00 

siny 

V 

y* 

00 

1.0 

00 

y" 

(3.64) 

z' 

-siny 

0.0 

cosy 

z” 

The  complete  transformation  is  given  by  the  product 


X 

+cosP 

— sinP 

0.0 

cosy 

00 

siny 

y 

= 

sinp 

cos(3 

0.0 

0.0 

1.0 

00 

z 

.  0.0 

0.0 

10. 

-siny 

00 

cosy 

cosy  cosfi 

-sinp 

siny  cosp 

X* 

= 

cosy  sinp 

cosp 

siny  sinp 

y" 

-siny 

00 
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z" 

This  transformation  matrix  is  used  by  CORFLD  to  determine  the 
global  radome  coordinates  from  the  local  antenna  coordinates  for  a 
gimbaled  antenna  rotated  by  (3  degrees  about  the  z-axis,  then  y 
degrees  about  the  new  y-axis.  With  this  computational  tool  the  fields 
that  are  received  by  an  antenna  within  the  core  of  the  radome  can  be 
determine  for  any  orientation  of  the  antenna  and  for  any  computed 
incident  field. 
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TV.  PROGRAM  EXECUTION  AND  DATA  FLOW 


A.  THE  RADOME  MODEL  DESIGN  PROCESS 

The  design  process  developed  by  this  thesis  is  represented  by 
Fig.  4.1.  The  radome  design  procedure  is  an  iterative  process  of 
selecting  radome  size  and  material  parameters,  radome  shaping 
functions,  and  placement  of  the  axis  origin  in  order  to  allow  for  an 
optimal  mesh  formation  by  EMRAD. 

The  flow  of  data  starts  when  RADSHP  is  used  to  design  the 
radome.  RADSHP  generates  a  material  file  which  is  ready  for  use  by 
EMRAD  and  a  structure  data  file  which  is  in  the  format  required  by 
CURVE  DIGITIZER.  CURVE  DIGITIZER  is  used  to  visually  check  the 
radome  structure  as  designed  in  RADSHP.  After  verifying  the 
structure  data  file,  the  program  EMCADIN,  developed  by  Connolly 
[Ref.  4],  is  used  to  translate  the  data  file  into  the  form  required  by 
EMRAD.  After  data  translation,  the  program  EMESH,  developed  by 
Morgan  [Ref.  3],  is  used  to  check  that  the  mesh  to  be  generated  in 
EMRAD  is  optimally  placed.  The  placement  of  the  origin  in  RADSHP 
when  generating  the  structure  data  file  can  be  modified  by  a 
placement  factor.  With  the  material  and  structure  files  are  prepared, 
EMRAD  is  executed  to  generate  the  interior  coefficients  of  expansion 
for  use  by  CORFLD.  A  sample  annotated  execution  through  the  data 
flow  follows. 
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a 


Figure  4.1  Radome  model  design  process 
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B.  RADSHP  EXECUTION 

RADSHP  is  a  menu-driven  program  that  accepts  data  from  the 
keyboard.  When  the  executable  version  of  RADSHP  is  available, 
RADSHP  is  invoked  by  typing 

>  RADSHP 


at  the  DOS  prompt.  The  program  will  load  itself,  then  provide 
introductory  information. 

WELCOME  TO  RADSHP 

THIS  PROGRAM  GENERATES  RADOME  MATERIAL  AND  SHAPE  DATA 
FILES  FOR  USE  WITH  EMRAD,  BASED  ON  USER  INPUT. 

RADSHP  prompts  for  data  file  names  for  the  material  and  structure 
output  files.  Sample  keyboard  inputs  are  provided.  User  inputs  are 
preceded  by  "==>”. 

TWO  FILES  ARE  REQUIRED  TO  HOLD  PROGRAM  OUTPUT  DATA. 

PLEASE  ENTER  THE  NAME  OF  THE  FILE  TO  BE  USED  TO 
HOLD  THE  PROGRAM  MATERIAL  PARAMETER  DATA.  THE  EXTENSION 
.MAT  WILL  BE  ADDED  AUTOMATICALLY,  I.  E.  YOURNAME.MAT 

==>  MATFLL 

PLEASE  ENTER  THE  NAME  OF  THE  FILE  TO  BE  USED 
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TO  HOLD  THE  PROGRAM  OUTPUT  DATA.  THE  EXTENSION  .DAT  WILL 
WILL  BE  ADDED  AUTOMATICALLY,  I.  E.  YOURNAME.DAT 


=>  STRFIL 

RADSHP  will  allow  the  construction  of  radomes  with  up  to  five  layers 
of  construction.  At  present,  each  of  these  layers  must  completely 
enclose  the  previous  layered  region. 

ENTER  THE  NUMBER  OF  LAYERS  (.LE.  5) 

==>  2 

RADSHP  prompts  for  the  parameters  of  each  layer  of  construction. 

Each  layer  of  the  radome  will  have  material  constants, 
an  outer  height  of  the  radome  layer,  H, 
a  outer  base  radius,  R,  and 

a  shape  function,  &  shape  factor  parameters  (if  appl.). 

The  first  "layer"  will  be  the  inner  core  of  the  radome, 
normally  filled  with  air. 

Layer  No:  1 

This  is  the  inner  core  of  the  radome 


Please  enter  the  radome  layer  height,  H, 
and  the  outer  base  radius,  R,  for  this  layer. 

==>  1.25,0.50 
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RADSHP  provides  a  menu  of  the  shaping  functions  available  for 
radome  construction.  The  curved  parabola  is  not  normally  used  for 
outer  layers  of  the  radome,  but  is  included  as  an  option  for  inner 
layer  construction. 

The  following  shaping  functions  are  available, 

1.  General  Ogive 

2.  Tangent  Ogive 

3.  von  Karman 

4.  Power  Series 

5.  Parabola,  curved 

6.  Parabola,  capped  by  pointed  apex 

Indicate  desired  shaping  function  by  entering  a  number 
1,2,  ...  or  6 


==>  2 


Each  shaping  function  requires  input  of  their  respective 
shaping  factors.  These  are  prompted  for  and  then  listed  for  each 
layer  as  the  radome  is  constructed.  The  tangent  ogive  and  the  von 
Karman  functions  have  shaping  factors  that  are  determined  by  the 
height  and  base  radius  of  the  radome  layer;  these  factors  are 
determined  and  listed.  The  factor  designation  follows  the  form 
adopted  in  the  equation  in  Chapter  II. 

The  tangent  ogive  shaping  function  determines  two 
factors  for  the  radome  construction.  These  are  calculated 
based  on  the  height  and  radius  of  the  radome  layer. 

These  factors  are  F4  =  6.50  and  F5  =  6.38 
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LAYER  SHAPE 

HEIGHT 

RADIUS 

FI 

F2 

F3 

F4 

F5 

F6 

F7 

1  TAN  Ogive 

1.25 

.50 

.00 

.00 

.00 

6.50 

6.38 

.00 

.00 

Layer  No:  2 

This  is  the  outer  layer  of  the  radome 


Please  enter  the  radome  layer  height,  H, 
and  the  outer  base  radius,  R,  for  this  layer. 

==>  1.30,0.55 

The  following  shaping  functions  are  available, 

1.  General  Ogive 

2.  Tangent  Ogive 

3.  von  Karman 

4.  Power  Series 

5.  Parabola,  curved 

6.  Parabola,  capped  by  pointed  apex 

Indicate  desired  shaping  function  by  entering  a  number 
1,2, ...  or  6 


==>  2 

The  tangent  ogive  shaping  function  determines  two 
factors  for  the  radome  construction.  These  are  calculated 
based  on  the  height  and  radius  of  the  radome  layer. 
These  factors  are  F4  =  6.42  and  F5  =  6.28 
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LAYER 

SHAPE 

HEIGHT 

RADIUS 

Fl 

F2 

F3 

F4 

F5 

F6 

F7 

1 

TAN  Ogive 

1.25 

.50 

.00 

.00 

.00 

6.50 

6.38 

.00 

.00 

2 

TAN  Ogive 

1.30 

.55 

.00 

.00 

.00 

6.42 

6.28 

.00 

.00 

After  all  the  structure  forming  data  is  entered,  RADSHP 
prompts  for  the  material  relative  permittivity  parameter  for  each 
layer  of  construction.  Layer  1  is  the  inner-most  region  and  the 
permittivity  will  normally  equal  unity. 

Layer  No:  1 

This  is  the  inner  core  of  the  radome 
Enter  Er  for  this  Layer: 

==>  1.0 

Layer  No:  2 

This  is  the  outer  layer  of  the  radome 
Enter  Er  for  this  Layer: 

==>  4.0 

With  the  structure  and  material  parameters  entered.  RADSHP 
prompts  for  headers  to  identify  the  radome  constructed.  This  header 
information  is  stored  in  the  material  data  file  and  will  follow  the  data 
files  throughout  the  execution  of  the  programs  in  this  package. 
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THE  HEADERS  ALLOW  THE  USER  TO  IDENTIFY  THIS 
SET  OF  DATA  FROM  ALL  OTHER  SETS. 

PLEASE  ENTER  HEADER  #1  (64  CHARACTERS  MAX) 


==>  This  is  a  test  data  set  for  RADSHP  example 


PLEASE  ENTER  HEADER  #2  (64  CHARACTERS  MAX) 


==>  Header  info  stays  with  data  set  throughout 


With  all  the  data  input  and  calculations  completed,  RADSHP 
reports  the  generation  of  the  material  and  structure  files  and 
terminates. 

Computations  and  output  are  completed,  the  data  file,  MA  I FlL. MAT 
holds  the  material  specifications  for  input  to  EMCAD. 

The  data  file,  STRFIL.DAT,  holds  the  radome  shape  data  for 
viewing  in  the  CAD  package,  CURVE  DIGITIZER. 

The  structure  file  required  for  input  to  EMCAD  can  be 
generated  by  the  executing  the  program  EMCADIN, 
using,  STRFIL.DAT,  as  input  to  EMCADIN, 

The  output  file  from  EMCADIN  will  be  in  proper  format 
for  input  to  EMCAD  for  the  structure  file. 

Stop  -  Program  terminated. 


> 
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C.  RADOME  CONSTRUCTION  VERIFICATION. 


The  construction  of  the  radome  by  RADSHP  is  verified  by 
entering  the  CURVE  DIGITIZER  package  and  displaying  the  structure 
data  file  generated  by  RADSHP.  Figure  4.2  is  a  copy  of  the  CURVE 
DIGITIZER  display  of  the  radome  shape  created  in  the  example 
execution  above.  Upon  review  and  acceptance  of  the  structure  data 
file,  EMCADIN  is  used  to  translate  the  data  file  into  the  form  required 
by  EMRAD.  EMCADIN  is  invoked  by  typing 

>EMCADIN 

at  the  system  prompt.  The  program  will  load  itself  and  displays  the 
greeting  heading  and  prompts  for  the  structure  data  file  as  input. 


WELCOME  TO  EMCADIN 

THIS  PROGRAM  ACCEPTS  PROPERLY  FORMATTED  INPUT  DATA 
FROM  CURVE  DIGITIZER  AND  CONVERTS  IT  TO  A  FORM 
WHICH  CAN  BE  USED  BY  EMCAD. 

PLEASE  PRESS  ENTER  TO  CONTINUE. 

INPUT  FILENAME 

INPUTFN  IS  THE  INPUT  FILENAME  OF  THE  FILE  CONTAINING 
THE  DATA  OF  INTEREST  WHICH  WILL  BE  INTERPOLATED  ON  AND 
CONVERTED  TO  A  FORM  WHICH  CAN  SUBSEQUENTLY  BE  USED  BY 
EMCAD 

PLEASE  ENTER  THE  NAME  OF  THE  INPUT  RLE.  THE  EXTENSION 
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Figure  4.2  Tangent  ogive  created  by  RADFLD 
as  displayed  by  CURVE  DIGITIZER 


OF  THE  FILENAME  MUST  BE  .DAT,  I.E.  FILENAME.DAT 


==>  STRFIL.DAT 


EMCADIN  then  prompts  the  user  for  the  data  file  name  to  hold  the 
structure  file  that  is  prepared  for  EMRAD  execution. 

OUTPUT  FILENAME 

OUTPUTFN  IS  THE  OUTPUT  FILENAME  OF  THE  FILE  CONTAINING 

THE  DATA  WHICH  HAS  BEEN  INTERPOLATED  AND  CONVERTED  TO  A 

FORM  WHICH  WILL  BE  USED  BY  EMCAD. 

PLEASE  ENTER  THE  NAME  OF  THE  OUTPUT  FILE.  THE  EXTENSION 

OF  THE  FILENAME  MUST  BE  INCLUDED  I.E.  FILENAME.DAT 

==>  STRFIL.STR 

There  are  several  data  files  used  within  these  packages,  the 
default  file  name  extensions  should  be:  ".MAT"  for  material  data  files, 
".DAT"  for  structure  data  files  suitable  for  CURVE  DIGITIZER,  and 
".STR"  for  structure  files  suitable  for  EMRAD. 

After  receiving  the  data  file  names  to  work  with,  EMCADIN 
prompts  for  angular  resolution  requirements.  EMCADIN  uses  linear 
interpolation  to  provide  the  capability  to  generate  output  data  with 
better  resolution  than  the  input  data.  The  default  resolution  of  data 
generation  by  RADSHP  is  one  thousand  points  in  the  x-z  plane,  using 
180  points  of  angular  resolution  in  EMCADIN  is  sufficient  for  good 
results. 
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RESOLUTION 


DELTHE  IS  THE  USER  INPUT  VALUE  OF  THE  DESIRED  THETA 
RESOLUTION  IN  DEGREES. 

PLEASE  ENTER  THE  DESIRED  DELTA  THETA  VALUE  IN  DEGREES. 

=>  180. 

Upon  completion  of  the  data  file  translation  by  EMCADIN,  the 
data  files  required  for  execution  by  EMRAD  are  complete  and  in  the 
proper  format.  MATFIL.MAT  is  the  example  material  file  and 
STRFIL.STR  is  the  example  structure  file.  Before  executing  EMRAD, 
the  efficiency  of  mesh  generation  to  be  expected  during  EMRAD  can 
be  previewed  by  using  the  program  EMESH  developed  by  Morgan. 
Figure  4.3  shows  the  mesh  generated  by  EMESH  for  the  example 
radome,  as  displayed  by  CURVE  DIGITIZER. 

D.  EMRAD  EXECUTION 

With  the  data  files  generated  and  an  optimal  mesh  solution,  the 
data  flow  is  ready  to  initiate  EMRAD.  EMRAD  is  a  menu-driven 
program  developed  by  Connolly  and  Morgan  and  has  been  slightly 
modified  to  provide  the  interior  spherical  coefficients  of  expansion 
for  use  by  CORFLD.  EMRAD  can  be  invoked  by  typing 

>EMCAD 

at  the  DOS  prompt.  The  program  will  load  itself  and  provide 
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Figure  4  3  Tangential  ogive  with  mesh  generated  by  EMESH 
as  displayed  by  CURVE  DIGITIZER 


introductory  information. 


******  WELCOME  TO  EMRAD  ****** 


THIS  PROGRAM  COMPUTES  EM  FIELD  SCATTERING  FROM  PENETRABLE 
BODIES  OF  REVOLUTION  USING  THE  COUPLED  AZIMUTHAL  POTENTIAL 
(CAP)  FORMULATION  IN  CONJUNCTION  WITH  A  VARIATIONAL  FINITE- 
ELEMENT  TECHNIQUE  AND  A  TRI-REGIONAL  UNIMOMENT  METHOD.  THE 
NUMERICAL  SOLUTION  IS  PERFORMED  USING  A  TWO-SWEEP  RICCATI 
TRANSFORM  WITH  CIRCUMFERENTIAL  MARCHING. 

UPDATED  VERSION  OF  ORIGINAL  U.C.  BERKELEY  CODE. 

MODIFICATION S  BY  M.A.  MORGAN  MAR  1987  -  DEC  1989 

E.M.  CONNOLLY  JUL  1987  -  APR  1988 
R.J.  VINCE  SEP  1989  -  DEC  1989  . 


PLEASE  PRESS  ANY  KEY  TO  CONTINUE. 


*******  emrad  INPUT  ******* 


EMRAD  ALLOWS  A  USER  TO  INPUT  THE  NECESSARY  INFORMATION 
ACCORDING  TO  THE  LEVEL  OF  HIS  EXPERTISE. 

1  NOVICE  LEVEL  -  GIVES  BRIEF  EXPLANATIONS/DESCRIPTIONS 

OF  INPUT  VALUES.  DESCRIBES  THE  REQUIRED 
FORMAT  FOR  THE  INPUT  VALUE.  GIVES 
TYPICAL  VALUES  WHERE  APPLICABLE. 

2  EXPERT  LEVEL  -  ASSUMES  THE  USER  IS  FAMILIAR  WITH  THE 

INPUT  PARAMETERS  AND  FORMATS.  SIMPLY 
PROMPTS  FOR  REQUIRED  INPUTS. 

3  EXIT  EMRAD 
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PLEASE  SELECT  THE  LEVEL  OF  EXPERTISE  BY  ENTERING  1,  2,  OR  3. 


==>  1 


The  names  of  the  data  files  generated  by  RADSHP  and  translated  by 
EMCADIN  are  entered  at  the  keyboard  when  prompted  by  EMRAD. 


EMRAD  REQUIRES  DATA  INPUT  FROM  TWO  FILES  PREPARED 
IN  ACCORDANCE  WITH  THE  USERS  MANUAL.  THESE  FILES  MUST 
CONTAIN  MATERIAL  PARAMETERS  AND  STRUCTURE  PARAMETERS, 
RESPECTIVELY. 

ENTER  MATERIAL  PARAMETER  FILE  (D:FTLENAME.EXTENSION)  : 


==>  MATFTL.MAT 


ENTER  STRUCTURE  DATA  FILE  (DrFILENAME.EXTENSION): 

==>  STRFIL.STR 

EMRAD  opens  both  files  and  displays  the  header  information  which 
RADSHP  placed  in  the  material  file. 

This  is  a  test  dataset  for  RADSHP  example 
Header  will  stay  with  dataset  throughout 

EMRAD  creates  several  datafiles.  The  first  will  be  the  diary  file  for 
the  execution  of  EMRAD.  This  file  will  have  the  ".OUT"  extension. 
The  second  data  file  generated  will  have  the  same  root  name  as  the 
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diary  file,  this  is  the  data  file  which  holds  the  spherical  coefficients  of 
expansion  for  the  core  field  generation.  This  second  file  will  have  the 
extension  ".INT"  and  will  be  used  as  input  to  CORFLD. 


THE  OUTPUT  DATA  FILE  IS  THE  FILE  CONTAINING  THE  OUTPUT 
RESULTS  OF  ALL  NUMERICAL  CALCULATIONS  CONDUCTED  BY  EMRAD. 
THE  FORMAT  FOR  THIS  INPUT  IS  FILENAME  ONLY.  NO  EXTENSION  IS 
REQUIRED  OR  DESIRED.  EMRAD  AUTOMATICALLY  APPENDS  AN 
EXTENSION  OF  .OUT  TO  YOUR  FILENAME. 

A  SECOND  OUTPUT  DATA  FILE  WILL  BE  CREATED  TO  STORE  THE 
COEFFICIENTS  FOR  THE  INTERIOR  CORE  EXPANSION.  EMRAD  WILL 
APPEND  AN  EXTENSION  OF  .INT  FOR  THIS  DATA  FILE. 

PLEASE  ENTER  THE  FILENAME  OF  THE  OUTPUT  DATA  FILES. 

==>  INTFIL 


EMRAD  was  developed  by  Morgan  and  Connolly  to  determine  the 
scattered  field  and  the  radar  cross  section  of  the  body  of  revolution 
under  study.  Four  output  data  files  are  created  for  input  to  separate 
graphing  routines. 


THE  OUTPUT  GRAPHICS  DATA  FILE  IS  THE  FILE  CONTAINING  THE 
OUTPUT  DATA  FROM  EMRAD  TO  BE  USED  AS  INPUT  TO  GRAPHING 
ROUTINES.  THE  FORMAT  FOR  THIS  INPUT  IS  FILENAME  ONLY.  NO 
EXTENSION  IS  REQUIRED  OR  DESIRED.  EMRAD  AUTOMATICALLY 
APPENDS  AN  EXTENSION  OF  .TMT,  .TMP,  .TET,  AND  .TEP  TO  YOUR 
FILENAME  AS  IT  PRODUCES  FOUR  OUTPUT  FILES  FOR  GRAPHICS. 

.TMT >  TM  INCIDENCE,  F-THETA 

.TMP >  TM  INCIDENCE,  F-PHI 


.TET >  TE  INCIDENCE,  F-THETA 

.TET >  TE  INCIDENCE,  F-PHI 


PLEASE  ENTER  THE  FILE  NAME  OF  THE  OUTPUT  DATA  FILE. 

=>  SCATGRF 

EMRAD  prompts  for  a  graphics  caption  to  be  included  with  the  data 
set,  this  caption  is  passed  along  to  CORFLD  and  is  included  with  each 
plotted  figure.  This  caption  affords  the  user  the  ability  to  identify 
this  data  set  from  all  other  sets. 

THE  GRAPHICS  CAPTION  IS  A  PERSONALIZED  CAPTION  ALLOWING  THE 
USER  TO  IDENTIFY  THIS  SET  OF  GRAPHS  FROM  ALL  OTHER  SETS. 

THE  MAXIMUM  LENGTH  OF  THIS  CAPTION  IS  64  CHARACTERS. 

NOTE:  WHEN  USED  WITH  THE  GRAPHICS  PACKAGE,  THE  PROGRAM  IS 
ABLE  TO  DIFFERENTIATE  BETWEEN  UPPER  CASE  AND  LOWER  CASE 
CHARACTERS. 

PLEASE  ENTER  ANY  GRAPHICS  CAPTION  YOU  DESIRE 
==>  Graphics  caption 

As  discussed  in  Appendix  A,  EMRAD  utilizes  a  semi-annular 
conformal  mesh  in  order  to  conduct  electromagnetic  field 
calculations.  EMRAD  prompts  the  user  for  a  minimum  and  maximum 
mesh  density.  EMRAD  conducts  the  analysis  of  the  radome  structure 
using  a  linear  variation  of  density  from  maximum  density  decreasing 
to  minimum  density. 

DMIN  AND  DMAX  ARE  PARAMETERS  OF  MESH  DENSITY  IN  TERMS  OF 
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ELEMENTS/INTERIOR  LAMBDA.  INPUT  VALUES  ARE  EXPECTED  TO  BE 
REAL,  I.E.  THE  DECIMAL  POINT  MUST  BE  INCLUDED. 

TYPICAL  VALUES  ARE  DMIN  =  10.  AND  DMAX  =  15. 


PLEASE  INPUT  DMIN 

=>  15 

PLEASE  INPUT  DMAX 

==>  18 

EMRAD  allows  for  analysis  of  up  to  five  incident  fields.  EMRAD 
prompts  for  the  number  of  incident  angles,  and  then  for  the  value  of 
each  angle,  in  degrees.  EMRAD  uses  the  poynting  vector  direction  for 
the  incident  angle.  This  translates  to  an  angle  of  180  degrees  for  an 
angle  incident  at  the  radome’s  appex,  135  degrees  for  an  angle 
incident  45  degrees  off  the  apex,  and  90  degrees  for  an  incident 
angle  off  the  beam  of  the  radome. 

THE  NUMBER  OF  INCIDENT  FIELD  ANGLES  IS  THE  TOTAL  NUMBER  OF 
INCIDENT  FIELDS  THAT  IMPINGE  ON  THE  OBJECT  OF  INTEREST. 

THIS  PROGRAM  ALLOWS  A  MAXIMUM  OF  FIVE  INCIDENT  FIELD 
ANGLES.  WHEN  ENTERING  YOUR  ANSWER  PLEASE  DO  NOT  INCLUDE  A 
DECIMAL  BECAUSE  THE  INPUT  MUST  BE  IN  INTEGER  FORMAT,  I.E. 


NA  =  3  ,  NA  =  1 


PLEASE  INPUT  THE  NUMBER  OF  INCIDENT  FIELD  ANGLES. 
==>  2 
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ENTER  INC  FLD  ANGLE  (DEG)  FOR  #  1 


=>  180 

ENTER  INC  FLD  ANGLE  (DEG)  FOR  #  2 
==>  135 

EMRAD  performs  scattering  calculations  by  starting  at  0=0,  and 
progressing  clockwise  to  0  =  180  in  equal  increments,  the  menu 
allows  the  user  to  determine  the  size  of  the  increments. 

THE  NUMBER  OF  SCATTERING  FIELD  THETA  POINTS  DETERMINES  THE 
SPACING  BETWEEN  THETA  POINTS  DURING  EMRAD  ITERATIONS  AND 
CALCULATIONS. 

DELTA  THETA  =  180  i  (NUMBER  THETA  POINTS  -  1)  SO... 

NUMBER  THETA  POINTS  =  37  . >  DELTA  THETA  =  5  DEGREES 

NUMBER  THETA  POINTS  =  19  - — >  DELTA  THETA  =  10  DEGREES 

WHEN  ENTERING  YOUR  ANSWER,  PLEASE  DO  NOT  INCLUDE  A  DECIMAL 
BECAUSE  THE  INPUT  MUST  BE  IN  INTEGER  FORMAT.  I.E.  NT  =  19 

PLEASE  INPUT  THE  NUMBER  OF  SCATTERING  FIELD  THETA  POINTS 

==>  19 

As  part  of  the  scattering  analysis,  EMRAD  allows  for  eight  angles  of 
user  observation. 
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THE  NUMBER  OF  PHI  ANGLES  IS  THE  TOTAL  NUMBER  OF  PHI  ANGLES. 
THIS  PROGRAM  ALLOWS  A  MAXIMUM  OF  THREE  PHI  ANGLES.  WHEN 
ENTERING  YOUR  ANSWER  PLEASE  DO  NOT  INCLUDE  A  DECIMAL 
BECAUSE  THE  INPUT  MUST  BE  IN  INTEGER  FORMAT.  I.E. 


NP  =  3  ,  NP  =  1 


PLEASE  INPUT  THE  NUMBER  OF  PHI  ANGLES. 

==>  0 

EMRAD  truncates  the  modal  analysis  of  the  electromagnetic  field 
infinite  fourier  expansions  for  the  incident  plane  waves  and  the 
interior  core  as  well  as  the  exterior  scattered  fields.  Mathematical 
estimates  of  truncations  are  provided  and  the  user  is  prompted  for 
truncation  limits.  Higher  truncations  may  increase  the  modelling 
accuracy. 

ENTER  MSTOP  (.LE.  13) 

ESTIMATED  ''MINIMUM"  VALUE  IS:  4 
==>  5 


Modal  truncations  for  the  number  of  spherical  harmonic  modes  in 
the  spherical  radome  core  and  external  scattering  region  are 
estimated.  The  user  is  prompted  for  truncation  limits  for  the 
spherical  harmonics. 
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ENTERING  NO.  OF  INTERNAL  AND  EXTERNAL  EXPANSION  MODES 


ESTIMATED  (KI*RMIN)  "MINIMUM"  INTERNAL  NI  =  4 

ENTER  N1  (.GE.  MSTOP): 


ESTIMATED  (KO*RMIN)  "MINIMUM"  EXTERNAL  N2  =  9 

ENTER  N2  (.GE.  MSTOP): 

==>  9 

Check  that  N1+N2  is  .LE.  35  ....  Otherwise  Abort 


==>  [RETURN] 


When  executing  the  numerical  solution,  EMRAD  requires  a  large 
amount  of  temporary  disk  storage  space.  The  amount  of  storage 
required  is  estimated  and  provided  to  the  user.  The  user  must  select 
a  disk  drive  that  has  adequate  memory  available  for  usage.  The  use 
of  a  ram  drive  configuration  can  increase  the  speed  of  execution  by 
EMRAD. 

*****  Temporary  Storage  of  Riccati  Matrices  ***** 

Estimated  Disk  Space  Needed  in  MegaBytes:  1.930000E+00 
ENTER  Hard  Drive  or  RAM  Disk  Letter  (  C,  D,  E  ... ): 

==>  E 
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With  all  of  the  data  and  parameters  specified,  EMRAD  begins  the 
numerical  analysis  of  the  electromagnetic  interaction  of  the  radome 
structure  and  the  specified  incident  plane  waves.  Execution  is 
monitored  by  reports  to  the  terminal  screen. 

INDEXING  PROGRAM  THROUGH  VALUES  OF  M 

M-LOOP  ....  M  =  0 

CALL  MESH 

ENTER  I-LOOP;  NO  STEPS  TO  COMPLETE:  38 
1=  1  NO  STEPS  TO  GO:  37 

CALL  LODER 
CALL  MARCH 

1=  2  NO  STEPS  TO  GO:  36 

EX  M-LOOP,  DATOUT 


***************  COMPLETED  ************** 


Stop  -  Program  terminated. 


Upon  termination  of  the  program  EMRAD,  the  coefficients  for  analysis 
of  the  radome  core  fields  have  been  stored  in  the  data  file,  ready  for 
CORFLD  execution. 


E.  CORFLD  EXECUTION 

CORFLD  is  a  menu  driven  program  that  accepts  data  from  the 
data  file  created  by  EMRAD  and  from  the  keyboard.  When  the 


62 


executable  version  of  CORFLD  is  available,  CORFLD  is  invoked  by 
typing 

>  CORFLD 

at  the  DOS  prompt.  The  program  will  load  itself,  then  provide 
introductory  information. 

WELCOME  TO  CORFLD 

CORFLD  ASSEMBLES  THE  FIELD  ON  THE  INTERIOR  OF  A  RADOME 
WITH  THE  COEFFICIENTS  SUPPLIED  BY  AN  EMRAD  OUTPUT  FILE. 


CORFLD  prompts  for  the  name  of  the  data  file  generated  by  EMRAD. 
Only  the  root  file  name  is  required,  CORFLD  attaches  the  ".INT"  to  the 
root  file  name. 


THE  OUTPUT  DATA  FILE  IS  THE  FILE  CONTAINING  THE  OUTPUT 
RESULTS  OF  INTERIOR  FIELD  CALCULATIONS  CONDUCTED  BY  EMRAD. 
THE  FORMAT  FOR  THIS  INPUT  IS  FILENAME  ONLY. 

NO  EXTENSION  IS  REQUIRED  OR  DESIRED. 

CORFLD  AUTOMATICALLY  APPENDS  AN  EXTENSION  OF  INT  TO  YOUR 
FILENAME. 

PLEASE  ENTER  THE  FILENAME  OF  THE  OUTPUT  DATA  FILE. 


==>  INTFIL 


Header  information  that  has  followed  the  data  flow  is  obtained  from 
the  ".INT"  data  file  by  CORFLD  and  displayed  on  the  screen.  Notice 
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that  the  graphics  caption  entered  during  EMRAD  execution  as  well  as 
the  headers  entered  during  RADSHP  execution  are  displayed. 


TEST  DATA  FOR  USER  EXAMPLE 

THIS  HEADING  FOLLOWS  THE  DATA  THROUGHOUT 

Graphics  caption 


CORFLD  reads  the  material  parameters  stored  in  the  ”.INT"  file  and 
displays  them  on  the  screen. 


SRE=  ( 1.000000, -5. 000000E-07) 

SRU=  ( 1.000000, -5. 000000E-07) 

KR  =  (1 .000000,- 1  .OOOOOOE-06) 

CORFLD  reads  the  number  of  incident  angles  evaluated  by  EMRAD 
and  displays  each  angle  as  the  poynting  vector  angle  of  incidence  and 
also  as  the  angle  of  incidence  to  the  radome  exterior.  The  radome 
incident  angle  is  displayed  so  that  the  user  can  use  the  angle  as  a 
reference  when  canting  the  planar  antenna. 

Incoming  Inc.  Angle  #  1 

EMRAD  Poynting  angle  of  incidence  =  180.000000 

Radome  incident  angle  =  O.OOOOOOE+OO 

Incoming  Inc.  Angle  #  2 

EMRAD  Poynting  angle  of  incidence  =  135.000000 

Radome  incident  angle  =  45.000000 

After  displaying  the  angles  of  incidence  available  for  consideration, 
EMRAD  prompts  the  user  for  the  angle  of  interest.  Only  one  angle 
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may  be  selected  for  each  execution  of  CORFLD.  In  order  to  consider 
other  angles,  CORFLD  must  be  re-executed  by  the  user. 

Select  Angle  of  Interest  By  Number 
=>  1 

The  angle  of  consideration  is  angle  #  1  =  O.OOOOOOE+OO  degrees 

After  an  angle  is  selected  for  analysis,  the  user  may  enable  a  rotation 
of  the  planar  array. 

The  planar  array  is  set  by  default  to  lie  in  the 
x-y  plane,  with  boresight  in  the  z-direction. 

Do  you  want  to  change  the  array  boresight  ? 

(  "Y"  or  "N"  ) 


==>  N 

CORFLD  reads  the  data  from  the  ".INT"  file  and  assembles  the  interior 
field  components. 

MSTART=  0  MSTOP=  5 
Reading  Coefficients  and  Constructing 
Fields  for  Each  Azimuthal  "m"  Mode 

Terminating  M-Loop  at  m=l  for  +/-  Z-Axis  Incidence 


Azimuthal  Mode  "m"  = 


0 


Reading  Coefficients 


Azimuthal  Mode  "m"  =  1 

Reading  Coefficients 
Generate  Fields  Across  Array 
SPHERICAL  FIELDS  COMPUTED 


When  the  interior  fields  have  been  computed,  the  user  may  check  for 
the  field  component  values  for  any  point  on  the  planar  array. 


Check  Values  of  Coordinates  and  E-Field  ?  (Y/N): 
==>  Y 

Enter  Array  Triple  Index  to  Check  (J,L,I) 

J=  x-point  1,2,..,  50 

L=  y-point  1,2,..,  50 

1=  1  for  TM,  2  for  TE 


XP,YP,ZP:  -5.115592E-02  -5.1 15592E-02  O.OOOOOOE+OO 
XC,YC,ZC:  -5.1 15592E-02  -5.1 15592E-02  O.OOOOOOE+OO 

IE-RADI  =  5.023 160E-01 
IE-THI  =  1.391270E-02 
IE-PHII  =  5.022836E-01 
IE-TOTI  =  7.104954E-01 


Look  at  Another  Point  ?  (Y/N): 
==>  N 


The  user  is  supplied  with  the  selection  menu  to  allow  access  to  any 
component  representation  of  the  interior  fields  across  the  planar 
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antenna  which  is  oriented  as  specified  by  the  user.  Component  phase 
and  magnitudes  of  the  interior  fields  as  well  as  the  total  field 
magnitudes  are  available. 


CORFU)  OUTPUT  GRAPHICS  PRESENTATION 
CORE  FIELD  REPRESENTATION  SELECTION  MENU 


Please  select  the  type  of  field  representation  of  interest 
CORE  will  display  the  following  field  representations: 


1.  TM  INCIDENT  FIELD,  E-THETA 

2.  TM  INCIDENT  FIELD,  E-PHI 

3.  TM  INCIDENT  FIELD,  E- RADIAL 

4.  TM  INCIDENT  FIELD,  E-TOTAL  FIELD  MAGNITUDE 

5.  TE  INCIDENT  HELD,  E-THETA 

6.  TE  INCIDENT  FIELD,  E-PHI 

7.  TE  INCIDENT  FIELD,  E-RADIAL 

8.  TE  INCIDENT  FIELD,  E-TOTAL  FIELD  MAGNITUDE 

9.  COMPARE  COMPUTED  FIELD  TO  INCIDENT  FIELD 

10.  CHANGE  ASPECT  RATIO 

0.  FINISHED  WITH  THIS  ANGLE  OF  INCIDENCE 
Indicate  your  selection  entering  "l","2”,...,"10",or  "0" 


==>  4 


When  the  user  has  made  a  selection  from  the  menu,  the  three 
dimensional  (3-D)  plotting  routine  is  called  and  the  user  is  allowed  to 
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adjust  the  view  angle  for  best  observation  of  the  3-D  plot,  see 
Fig.  4.4, 


CALL  PLOTTING  ROUTINE 

A  black  and  white  'stick'  representation  as  well  as  a  color  3-D  fill 
representation  are  available  to  the  user.  The  ’stick’  representation 
generates  the  clearest  output  for  printing  and  is  used  for  all  the  data 
presented  here. 

Enter  0  for  3D-stick  or  1  for  3D-fiil  option 
==>  0 

(**  SEE  FIGURE  4.4  **) 

Default  view  angles  are:  phi=-45.,theta=70. 

ENTER  NEW  ANGLES  PHI  &  THETA  (DEG)  OR  (9,9)  WHEN  DONE 
==>  9,9 

When  the  user  is  finished  with  an  electric  field  component,  the  new 
view  angle  "9,9"  is  entered  and  the  user  is  returned  to  the  core  field 
selection  menu. 

CORE  FIELD  REPRESENTATION  SELECTION  MENU 


9.  COMPARE  COMPUTED  FIELD  TO  INCIDENT  FIELD 
Indicate  your  selection  entering  "1  ",”2",.. .,"10", or  "0" 
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The  datafile  used  to  generate  fields  is  intfil.INT 
TEST  DATA  FOR  USER  EXAMPLE 
THIS  HEADING  FOLLOWS  THE  DATA  THROUGHOUT 
Graphics  caption 

TM  INCIDENT  FIELD,  E-TOTAL  FIELD  MAGNITUDE 


win.  Z  =  +.631;  Max.  Z  =  +.BZ5;  average  Z  =  +.656 
Inc.  A  =  +.000;  Plane  G  =  +.000;  View:  phi=+45.0  theta=  +  75.0 


Figure  4  4  3-D  plot  of  E-total  field  magnitude  for  example  ogive  radome 
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When  the  user  wishes  to  see  the  amount  of  perturbation  caused  to 
the  incident  plane  wave  by  the  radome,  selection  "9"  is  entered, 


=>  9 

Selection  "9"  sends  the  user  into  Menu  2,  for  core  field  comparison  to 
the  ideal  theoretical  plane  wave  which  would  be  incident  on  the 
antenna  without  the  radome  present. 


CORE  FIELD  REPRESENTATION  SELECTION  MENU  2 


Please  select  the  type  of  field  representation  of  interest 
CORE  will  display  the  following  field  representations: 

1 .  TM  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field 

2.  TM  INCIDENT  FIELD,  Scaled  Error  component 

3.  TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field 

4.  TE  INCIDENT  FIELD,  Scaled  Error  component 


Indicate  your  selection  entering  ’T","2","3",  or  "4’' 
==>  2 

CALL  PLOTTING  ROUTINE 

Enter  0  for  3D-stick  or  1  for  3D-fiH  option 

==>  0 

Default  view  angles  are:  phi=-45.,theta=70. 

(**  SEE  FIGURE  4.5  **) 


70 


Figure  4.5  3-D  plot  of  scaled  error  component  for  example  ogive  radome 


ENTER  NEW  ANGLES  PHI  &  THETA  (DEG)  OR  (9,9)  WHEN  DONE 


==>  9,9 

CORFLD  OUTPUT  GRAPHICS  PRESENTATION 
CORE  FIELD  REPRESENTATION  SELECTION  MENU 

0.  FINISHED  WITH  THIS  ANGLE  OF  INCIDENCE 
Indicate  your  selection  entering  "1", "2" ,...,"10", or  "0" 

==>0 

*********  CORE  ANALYSIS  COMPLETED  ************** 

Stop  -  Program  terminated. 

When  the  user  is  finished  with  the  angle  of  incidence  selected,  "0"  is 
selected  at  the  selection  menu  and  the  program  execution  is 
terminated. 
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V.  PROGRAM  VALIDATION  AND  RESULTS 


In  order  to  validate  the  execution  of  EMRAD  and  CORFLD,  the 
tangent  ogive  radome  generated  in  Chapter  IV  was  modified  to  have 
a  relative  permittivity  equal  to  that  of  air,  er=1.0.  The  field 
coefficients  were  generated  by  EMRAD  and  the  core  fields  were 
computed  by  CORFLD  and  compared  to  the  unperturbed  incident  field 
as  listed  in  Chapter  III. 

Figures  5.1  and  5.2  show  the  magnitude  and  phase  of  the 
scaled  dot  product,  as  defined  in  Chapter  III,  between  the  computed 
and  theoretical  values  for  the  TM  and  TE  incident  fields,  respectively, 
for  a  "nose-on"  incident  angle.  Each  comparison  shows  a  94.9% 
magnitude  matching  with  an  almost  constant  2.59°  phase  differential 
across  the  array.  Figure  5.3  exhibits  the  scaled  error  component  for 
the  TM  and  TE  incident  fields,  with  an  average  magnitude  of  error 
component  of  0.4%.  The  comparisons  for  this  "nose-on"  angle  of 
incidence  show  that  there  is  good  agreement  between  the  theoretical 
and  calculated  field  components.  It  is  observed  that  the  TM  and  TE 
components  displayed  are  90°  out  of  phase  with  each  other.  This  is 
expected  for  "nose-on"  incidence  since  the  TM  wave  has  its  direction 
component  in  the  -x  direction,  while  the  TE  wave’s  direction 
component  is  in  the  y  direction. 

Figure  5.4  shows  the  magnitude  and  phase  of  the  scaled  dot 
product  between  the  computed  and  theoretical  values  for  the  TM 
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Tlie  datafile  used  to  generate  fields  is  airfil.INT 
Shaping  function  =  TAN  ogive 
Naterial  of  formation  is  air*  Er  =  1 


Ttl  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  NAGN. 


min.  Z  =  +.982;  max.  Z  =  +.99S;  average  Z  =  +.949 
Inc.  A  =  +.000;  Plane  G  =  +.000;  Ulew:  pl»i=+45.0  theta=+75.0 


The  datafile  used  to  generate  fields  is  airfll.IHT 
Shaping  function  =  TAN  ogive 
Naterial  of  formation  Is  air,  Er  =  1 


TN  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  PHASE 


Inc.  A  =  +.000;  Plane  G  =  +.000;  Ulew:  plii=  +  80. 0  tl»eta=  +  7G.0 


Figure  5.1  Dot  product  magnitude  and  phase  for  TM 
computed  fields  determined  for  radome  composed  of  air 
compared  with  theoretical  plane  waves 
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The  datafile  used  to  generate  fields  is  airfil.INT 
Shaping  function  =  TAN  ogive 
Material  of  fornation  is  air,  Er  =  1 


TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  MAGN. 


min.  Z  -  +.982;  max.  Z  =  +.996;  average  Z  =  +.949 
Inc.  A  =  +.000;  Plane  G  =  +.000;  Uiewl  phi=  +  45.0  theta=+75.0 


The  datafile  used  to  generate  fields  is  airfil.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  air,  Er  =  1 

TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  PHASE 


Figure  5  2  Dot  product  magnitude  and  phase  for  TE 
computed  fields  determined  for  radome  composed  of  air 
compared  with  theoretical  plane  waves 
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The  datafile  used  to  generate  fields  is  airfil.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  ait,  Er  =  1 

TE  INCIDENT  FIELD,  Scaled  Error  component 
0. 10  r- 


0.  05 


-0.05 


min.  Z  =  +.000;  wax.  Z  =  +.023;  average  Z  =  +.004 
Inc.  A  =  +.000; Plane  G  =  +.000;  View:  phi=+10.0  theta=+80. 0 


The  datafile  used  to  generate  fields  is  airfil.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  air,  Er  =  1 

TM  INCIDENT  FIELD,  Scaled  Error  component 

0. 10  r- 


0.05 


min.  Z  =  +.000;  max.  2  =  +.023;  average  Z  =  +.004 
Inc.  A  =  +.000; Plane  G  =  +.000;  View:  phi=+75.0  theta=+80. 0 


Figure  5.3  Scaled  error  component  for  TM  and  TE 
computed  fields  determined  for  ra dome  composed  of  air 
compared  with  theoretical  plane  waves 
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The  datafile  used  to  generate  fields  is  airfil.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  air,  Er  =  1 

TM  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  MAGN. 


win.  Z  -  *.  540;  max.  Z  =  +1.66;  average  Z  =  +.874 
Inc.  A  =  +60.  0;  Plane  G  =  *60.0;  Uiew:  phi=  +  45.0  tlieta=  +  75.0 

The  datafile  used  to  generate  fields  is  airfil.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  air,  Er  =  1 


TM  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  PHASE 


min.  Z  =  -27.9;  max.  Z  ~  +60. i;  average  Z  =  -2.84 
Inc.  A  =  +60.0;  Plane  G  =  +60.0;  Uieu:  phi=  +  45.0  lhela=  +  75.0 


Figure  5 A  Dot  product  magnitude  and  phase  for  TM  computed  fields 
determined  for  radome  composed  of  air  compared  with  theoretical 
plane  waves  for  incident  angle  of  60  degrees 
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incident  field,  for  an  incident  angle  60°  from  the  "nose-on"  aspect. 
The  magnitude  of  the  scaled  dot  product  between  the  computed  and 
theoretical  TM  incident  wave  shows  modal  variation  around  the 
antenna  plane.  The  magnitude  of  the  computed  value  in  the 
direction  of  the  theoretical  plane  wave  varies  between  54.0%  and 
166%,  with  an  average  magnitude  of  87.4%.  The  phase  component  of 
the  dot  product  varies  across  the  plane  from  -27.9°  to  +60.1°,  with 
an  average  phase  of  -2.84°.  The  maximum  phase  differential  occurs 
at  the  edge  of  the  antenna  plane  closest  to  the  point  of  wave 
incidence,  as  the  incident  wave  is  defined  as  arriving  in  the  -x  -  z 
axis.  This  variation  seen  in  the  magnitude  and  phase  of  the  scaled 
dot  product  in  the  direction  of  the  theoretical,  undisturbed  plane 
wave  is  confirmed  in  the  scaled  error  component  of  Fig.  5.5  with  an 
amplitude  variation  from  8.4%  to  88.6%,  averaging  at  23.6%  of  the 
computed  field  evaluated  as  an  error  term  perpendicular  to  the 
theoretical  fields  expected. 

Figures  5.6  and  5.7  exhibit  the  behavior  of  the  calculated  TE 
60°  incident  plane  wave,  which  is  similar,  but  less  in  magnitude 
compared  to  the  calculated  TM  wave.  Less  modal  variation  is  seen 
along  the  exhibited  plane  of  the  antenna.  The  scaled  error 
component  has  a  symmetry  about  the  x-axis.  The  scaled  error 
amplitude  varies  form  0.3%  to  33.1%,  with  an  average  of  12.5%. 

The  form  of  the  error  components  as  analyzed  for  this  radome 
constructed  of  air  identifies  a  phase  error  in  the  determination  of  the 
computed  fields.  The  data  computations  of  the  EMRAD  subroutines 
have  been  previously  validated  by  Morgan.  These  validations 
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The  datafile  used  to  generate  fields  is  airfil.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  air,  Er  =  1 


TM  INCIDENT  FIELD,  Scaled  Error  component 


min.  Z  =  +.884;  max.  Z  =  +.886;  average  Z  =  +.Z36 
Inc.  A  =  +60. 8;  Plane  G  =  +60.0;  View:  phi=+45.0  theta=  +  75.  0 


Figure  5  5  Scaled  error  component  fo.  TM  computed  fields 
determined  for  radome  composed  of  air  compared  with 
theoretical  plane  waves  for  incident  angle  of  60  degrees 


The  datafile  used  to  generate  fields  is  airfil.INT 
Shaping  function  =  TAN  ogive 
Naterial  of  formation  is  air,  Er  =  1 

TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  MAGN. 

2.00  tr 


1.50 


min.  Z  =  +.  860;  max.  Z  =  41.34;  average  Z  =  4.981 
Inc.  A  =  460.0;  Plane  G  =  460.0;  Uiew:  phi=445.8  theta=475.0 


The  datafile  used  to  generate  fields  is  airfil.INT 
Shaping  function  =  TAN  ogive 
Material  of  formation  is  air,  Er  =  1 

TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  PHASE 

Z7.00  E- 


min.  Z  =  -18.2;  max.  Z  =  44.96;  average  Z  =  -.438 
Inc.  A  =  460. 0;  Plane  G  =  460.0;  Uiew:  phi  =  +  45.0  theta=475.0 

Figure  5.6  Dot  product  magnitude  and  phase  for  TE  computed  fields 
determined  for  radome  composed  of  air  compared  with  theoretical 
plane  waves  for  incident  angle  of  60  degrees 
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Die  datafile  used  to  generate  fields  Is  airfil.INI 
Shaping  function  =  TAN  ogive 
Material  of  fornation  Is  air,  Er  =  1 


TE  INCIDENT  FIELD, 
0.50  r~ 


Scaled  Error  coMponent 


00 


win.  Z  =  +.003;  Max.  Z  =  +.331;  average  Z  =  +.125 
Inc.  A  =  *80.0;  Plane  G  =  +60.  0,'  Uieu:  phl=  +  20. 0  tlieta=  +  75.0 


Figure  5  7  Scaled  error  component  for  TM  computed  field 
determined  for  radome  composed  of  air  compared  with 
theoretical  plane  waves  for  incident  angle  of  60  degrees 
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considered  the  magnitude  and  phase  of  the  scattered  fields  [Refs.  3 
and  4],  and  the  magnitude  of  the  core  fields  [Ref.  10]  along  the  main 
coordinate  axes.  This  work  is  the  first  to  plot  both  the  phase  and 
magnitude  variations  across  the  entire  two-dimensional  plane  and 
identifies  the  possible  need  for  increased  computational  accuracy  in 
the  numerical  routines  of  EMRAD.  There  is  also  the  possibility  of 
unresolved  programming  errors  in  CORFLD  which  should  be  further 
investigated. 

B.  COMPUTATION  RESULTS  FOR  VARIOUS  RADOME  MATERIALS 

EMRAD  and  CORFLD  executions  for  the  tangent  ogive  standard 
established  in  Chapter  IV  were  run  for  varying  material  permittivity 
according  to  the  material  properties  listed  in  Table  2.1.  The 
comparisons  of  the  computed  fields  to  the  unperturbed  theoretical 
incident  plane  waves  are  included  only  for  the  "nose-on"  incident 
angle  which  has  good  results  for  the  ideal  "air”  construction. 
Fiber-reinforced  glass,  with  er=4.0;  glass  ceramic,  with  er=5.5;  and 
refractory  oxide,  with  er=8.0,  are  included.  The  form  of  field 
comparisons  are  similar  in  shape  for  each  radome  considered,  only 
the  magnitude  of  the  plots  varies.  As  expected  from  the  discussion 
in  Chapter  II,  the  attenuation  of  the  incident  plane  wave,  received  at 
the  plane  of  the  antenna  in  the  core,  increases  with  increasing 
relative  permittivity.  When  the  radome  was  constructed  of  a;r,  the 
computed  magnitude  of  each  incident  wave  was  94.9%.  When 

fiber-reinforced  glass  construction  is  considered.  Figures  5.8,  5.9,  and 
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5.10,  the  average  magnitude  of  the  computed  wave  is  60.1%  and  the 
average  scaled  error  component  is  22.7%.  For  the  radome 
constructed  of  glass  ceramic.  Figures  5.11,  5.12,  and  5.13,  the  average 
magnitude  is  66.9%  and  the  average  scaled  error  component  is  13.3%. 
For  the  refractory  oxide  radome.  Figures  5.14,  5.15,  and  5.16,  the 
average  magnitude  of  received  field  is  87.6%  and  the  average  scaled 
error  component  is  20.6%. 


83 


The  datafile  used  to  generate  fields  is  P4FIL. INI 
Shaping  function  =  TAN  ogive 

Naterial  of  formation  is  fiber-reinforced  glass 


TN  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  MAGN. 


win.  Z  =  +.475;  max.  Z  =  +.811;  average  Z  =  +.601 
Inc.  A  =  +.000;  Plane  C  =  +.000;  Uiew:  phi=+35.0  theta=+70. 0 


The  datafile  used  to  generate  fields  is  p4fil.INT 
Shaping  function  =  TAN  ogive 

Material  of  formation  is  fiber-reinforced  glass 


T n  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  PHASE 


ormin.  Z  =  -71.1;  max.  Z  =  -27.4;  average  Z  =  -51.8 
Inc.  A  =  +.000; Plane  G  =  +.000;  View:  phl  =  +  10. 0  theta=  +  85.0 


Figure  5.8  Dot  product  magnitude  and  phase  for  TM  computed  fields 
determined  for  radome  composed  of  fiber-reinforced  glass  compared 
with  theoretical  p!?ne  waves  for  incident  angle  of  0  degrees 
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The  datafile  used  to  generate  fields  is  p4fil.INT 
Shaping  function  =  TAN  ogive 

Naterial  of  forMation  is  fiber-reinforced  glass 
T M  INCIDENT  FIELD,  Scaled  Error  component 


min.  Z  =  +.813;  max.  Z  =  +.436;  average  Z  =  *.22? 
Inc.  A  =  ♦  . 080; Plane  G  =  +.800;  Uleu:  phi =+95.0  theta=+70. 0 


The  datafile  used  to  generate  fields  is  p4fiI.INI 
Shaping  function  =  TAN  ogive 

Naterial  of  formation  is  f iber-rei nforced  glass 
TE  INCIDENT  FIELD,  Scaled  Error  component 


min.  Z  =  +.813;  max.  Z  =  +.436;  average  Z  =  +.227 
Inc.  A  =  +.000;  Plane  G  =  +.000:  View:  pl,i=  +  10.0  tl»eta=470.  0 

Figure  510  Scaled  error  component  for  TM  and  TE  computed  fields 
determined  for  radome  composed  of  fiber-reinforced  glass  compared 
with  theoretical  plane  waves  for  incident  angle  of  0  degrees 
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The  datafile  used  to  generate  fields  is  c55f 1 1  INT 
TEST  DATA  FOR  THESIS 


Radome  is  tangential  ogive, 
constructed  of  ceranic,  Er  =  5.5 

TM  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  NAGN. 


The  datafile  used  to  generate  fields  is  c55fil.INT 
TEST  DATA  FOR  THESIS 
Radone  is  tangential  ogive, 
constructed  of  ceramic,  Er  =  5.5 

TM  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  PHASE 


min.  Z  =  -77.5;  max.  Z  =  -47.5;  average  Z  =  -61.9 
Inc.  A  =  ♦ .000;Plane  G  =  ♦.000;  Uieu:  phi=+45.0  theta=^75.0 


Figure  5  1  1  Dot  product  magnitude  and  phase  for  TM  computed  fields 
determined  for  radome  composed  of  glass  ceramic  compared 
with  theoretical  plane  waves  for  incident  angle  of  0  degrees 
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The  datafile  used  to  generate  fields  is  c55fil.INT 
TEST  DATA  FOR  THESIS 
Radome  is  tangential  ogive, 
constructed  of  ceram icj  Er  =  5.5 

TE  INCIDENT  FIELD,  Computed  Field  dot  Idea)  Field,  MAGN. 


min.  Z  =  *.668;  max.  Z  =  +.883;  average  Z  =  +.669 
Inc.  A  =  +.000;Plane  G  =  +.000;  Uieu:  phl=+45.0  theta=+75.0 


Tlie  datafile  used  to  generate  fields  is  cSSfil.INT 
TEST  DATA  FOR  THESIS 
Radome  is  tangential  ogive, 
constructed  of  ceramic,  Er  =5.5 

TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  PHASE 


min.  Z  =  -77.5."  max.  Z  =  -47.5;  average  Z  =  -61.9 
Inc.  A  =  +.000;Plane  G  =  +.008;  Uieu:  pl»l=+45.8  tlieta  =  +75.0 


Figure  512  Dot  product  magnitude  and  phase  for  TE  computed  fields 
determined  for  radome  composed  of  glass  ceramic  compared 
with  theoretical  plane  waves  for  incident  angle  of  0  degrees 
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The  datafile  used  to  generate  fields  is  c55fil.INT 
TEST  DATA  FOR  THESIS 
Radone  is  tangential  ogive, 
constructed  of  ceramic.  Er  =5.5 

TN  INCIDENT  FIELD.  Scaled  Error  component 


mii<.  Z  =  *.005;  max.  Z  =  ♦.312;  average  Z  =  ♦ .  133 
Inc.  A  =  ♦ .000; Plane  G  =  *.000;  Uieu:  phi =*80.0  theta =+70.0 


The  datafile  used  to  generate  fields  is  c55fil.INT 
TEST  DATA  FOR  THESIS 
Radone  is  tangential  ogive, 
constructed  of  ceramic.  Er  =  5.5 

TE  INCIDENT  FIELD.  Scaled  Error  component 
0.  50  r- 


mln.  Z  =  *.0051  max.  Z  =  ♦,31Z»'  average  Z  =  +.133 
Inc.  A  =  +.000;  Plane  G  =  t.000;  Mieu:  phi=  +  10. 0  theta=  +  70. 0 


Figure  513  Scaled  error  component  for  TM  and  TE  computed  fields 
determined  for  radome  composed  of  glass  ceramic  compared 
with  theoretical  plane  waves  for  incident  angle  of  0  degrees 
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The  datafile  used  to  generate  fields  is  0X8FIL.INT 
Shaping  function  =  TAN  ogive 

Material  of  formation  is  refractory  oxide,  Er  =  8 

TM  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  MAGN. 

2.80  t- 

1.S0  I- 


1.88 


min.  Z  =  *.6341  max.  Z  =  *1.081  average  Z  =  *.876 
Inc.  A  =  ♦  .800  IP  lane  G  =  *.000:  Uieut  phi=*4S.0  theta=*75.0 


Tlte  datafile  used  to  generate  fields  Is  0X8FIL.INT 
Shaping  function  =  TAN  ogive 

Materia]  of  formation  is  refractory  oxide,  Er  =  8 


TM  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  PHASE 


Figure  5-M  Dot  product  magnitude  and  phase  for  TM  computed  fields 
determined  for  radome  composed  of  refractory  oxide  compared 
with  theoretical  plane  waves  for  incident  angle  of  0  degrees 
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The  datafile  used  to  generate  fields  is  0X8FIL.INT 
Shaping  function  =  TAN  ogive 

Material  of  formation  is  refractory  oxidet  Er  =  8 

TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  NAGN. 

2.00  |- 

1.50  §- 


1.00 


min.  Z  =  ♦ .634;  max.  Z  =  +1.08;  average  Z  =  +.876 
Inc.  A  =  *.000;Plane  G  =  +.000;  Uieu:  phi=*5S.0  theta=+75.0 


The  datafile  used  to  generate  fields  is  0X8FIL.INT 
Shaping  function  =  TAN  ogive 

Hater ia I  of  formation  is  refractory  oxide,  Er  =  8 


TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  PHASE 


Inc.  A  =  +.000;Plane  G  =  +.000;  Uieu:  phi=*85.0  theta=-»80.0 


Figure  5.15  Dot  product  magnitude  and  phase  for  TE  computed  fields 
determined  for  radome  composed  of  refractory  oxide  compared 
with  theoretical  plane  waves  for  incident  angle  of  0  degrees 
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The  datafile  used  to  generate  fields  is  0X8FIL.1NT 
Shaping  function  =  TAN  ogive 

Material  of  formation  is  refractory  oxide,  Er  =  8 


TM  INCIDENT  FIELD,  Scaled  Error  component 


min.  Z  =  *.006;  max.  Z  =  +.560;  average  Z  =  +.206 
Inc.  A  =  ♦.000iPIane  G  =  ♦.000;  Uieui  phi=+85.0  theta=+75.0 


The  datafile  used  to  generate  fields  is  OXBFIL.INT 
Shaping  function  =  TAN  ogive 

Material  of  formation  is  refractory  oxide,  Er  =  8 


TE  INCIDENT  FIELD,  Scaled  Error  component 


min.  Z  =  ♦.006;  max.  Z  =  +.560;  average  Z  =  *.206 
Inc.  A  =  ♦.800;Plane  G  =  ♦.000;  Uieu:  phi=*25.0  theta=*7S.0 


Figure  5.16  Scaled  error  component  for  TM  and  TE  computed  fields 
determined  for  radome  composed  of  refractory  oxide  compared 
with  theoretical  plane  waves  for  incident  angle  of  0  degrees 
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VI.  CONCLUSIONS 


This  thesis  has  presented  a  useful  set  of  programs  to 
numerically  model  a  radome,  designed  of  specified  shape, 
construction,  and  material,  and  to  evaluate  the  perturbation  of  the 
electromagnetic  field  as  it  passes  through  the  modelled  radome  to 
the  antenna  within  the  radome  structure.  The  radome  design 
methodology  was  discussed  in  Chapter  II,  along  with  the  radome 
shaping  functions,  wall  construction,  and  choice  of  materials  for 
radome  construction.  The  program  RADSHP  is  introduced  and  the 
shaping  function,  as  well  as  the  shaping  factors  for  each  function  are 
explained.  Chapter  III  developed  the  electromagnetic  theory 
required  to  assemble  the  field  components  within  the  interior 
radome  core,  using  the  coefficients  provided  by  the  EMRAD  program. 
The  theory  required  for  transforming  the  coordinates  of  the  canted 
antenna  plane  to  the  core  coordinates  for  an  assumed  two-axis 
gimbal  was  also  discussed.  Chapter  IV  explained  the  data  flow  for 
the  radome  model  design  process  by  using  an  example  of  the 
program  executions  for  the  design  of  a  monolithic  tangent  ogive,  with 
computational  analysis  by  EMRAD  and  field  component  generation  by 
CORFLD.  The  verification  of  the  radome  design  and  optimal  placing  of 
the  core  origin  by  the  use  of  CURVE  DIGITIZER  software  is  also 
discussed.  Chapter  V  discussed  the  computational  errors  discovered 
in  the  3-D  plotting  of  the  field  components.  The  scaled  magnitude  of 
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the  computed  perturbed  plane  waves  for  the  "nose-on"  incident 

angle  agreed  within  95%  of  the  theoretical  value  and  the  phase  of  the 
computed  field  was  an  almost  constant  2.59°  out  of  phase  with  the 
theoretical  incident  plane  wave.  The  errors  increase  when  the  angle 
of  incidence  is  off  the  "nose-on"  aspect,  with  modal  amplitude 

variations  from  54%  to  166%  of  the  theoretical  plane  wave.  Further 
diagnostics  are  required  to  determine  the  source  of  this  error. 

This  is  only  the  first  stage  in  the  development  of  a 

mathematical  model  for  the  interaction  between  the  radome  and  the 
antenna  system  that  the  radome  is  designed  to  enclose.  The  model, 
in  its  present  form,  does  not  consider  the  presence  of  the  antenna 

structure  within  the  radome  core.  Further  work  is  needed  to 
determine  the  effect  of  the  antenna  structure  and  its  interaction  with 
the  radome.  When  such  a  model  is  developed  the  predicted  data  can 
be  compared  to  data  collected  by  a  physical  antenna  system.  The 
culmination  of  this  effort  could  be  the  analytic  development  of 
corrective  lenses,  as  referred  to  in  Chapter  II,  to  correct  or  minimize 
the  perturbation  of  the  incident  electromagnetic  field  caused  by  the 
radome  structure. 
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Appendix  A 


EMRAD  FIELD  THEORY 

The  program  EMRAD  uses  the  unimoment  method  mentioned  in 
Chapter  II,  as  developed  by  Morgan  and  Mei  [Ref.  5]  to  decouple  the 
interior  and  exterior  spherical  solutions.  The  expansion  coefficients 
for  the  electromagnetic  fields  within  the  radome  structure  are 
generated  as  part  of  the  scattering  solution.  CORFLD  uses  these 
coefficients  for  construction  of  the  interior  fields.  The  theory 
required  by  CORFLD  to  formulate  the  electromagnetic  fields  within 
the  radome  is  discussed  in  Chapter  II. 

The  code  developed  by  Morgan  for  EMRAD  uses  an  optimized 
variational  finite-element  algorithm,  in  conjunction  with  a 
tri-regional  unimoment  method,  to  provide  scattering  solutions  for 
each  of  multiple  incident  fields  impinging  on  the  specified  structure. 
EMRAD  represent  a  significant  improvement  in  computational  speed 
and  versatility  over  computer  codes  which  find  the  solution  to  the 
electromagnetic  scattering  and  penetration  problem  by  use  of  surface 
current  integral  equation  formulations.  This  advantage  results 
because  of  the  highly  sparse  matrices  embodied  in  the  finite  element 
solution  vis-a-vis  full  matrices  found  using  an  integral  equation 
formulation.  An  overview  of  the  finite  element  scattering  algorithm, 
along  with  the  concept  of  coupled  azimuthal  potentials  and  the 
unimoment  method  are  included  in  this  appendix.  Further 
information  is  available  in  references  3,  4,  5,  6,  and  7. 
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A.  COUPLED  AZIMUTHAL  POTENTIALS 

In  the  radome  construction  the  material  constitutive 
parameters,  e(r)  and  p(r),  are  invariant  to  the  azimuthal  coordinate 
due  to  the  axial  symmetry  of  the  radome.  This  axisymmetry  allows 
implementation  of  vector  field  generation  via  coupled  azimuthal 
potentials,  as  developed  by  Morgan,  Chang  and  Mei  [Ref.  7]. 

For  ease  of  solution  and  representation,  a  normalized 
cylindrical  coordinate  system  is  adopted,  as  defined  by  (R  ,  Z  ,  <j>)  = 
(k0  p,  kQ  z,  4>)  where  (p,  z,  <}>)  are  standard  circular  cylindrical 
coordinates  and  k0  =  2jiA0  is  the  free-space  wavenumber  of  the 
time-harmonic  field.  The  equations  in  this  form  are  related  to  the 
standard  spherical  coordinates  (r,  6,  <)>)  as  shown  in  Fig.  A.l  using  the 
substitutions 


R  =  kG  r  sin0 


(A.l) 


and 


Z  =  k0  r  cos0  .  (A. 2) 

The  total  electromagnetic  fields  are  decomposed  into  azimuthal 
modes  by  an  exponential  Fourier  series  in  the  <J> -coordinate  such  that 


E(R,Z,«  =  £  em(R.Z)  ei"1* 


(A.3) 


m= 


and 


n„H(R,Z,<»  =  £  hm(R,Z)  eJ1"0  . 


(A. 4) 


m=- 
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Figure  A.l  Rectangular  (x,  y,  z)  spherical  (r;e,  40,  and 
normalized  cylindrical  (R,  Z,8)  coordinate  systems 
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where  r\0  =  120  n  £3  is  used  to  normalize  the  magnetic  fields,  so  that 

h  =  t\0  H. 

Using  the  axial  symmetry  of  the  radome  construction,  the 
modal  fields  decouple  when  the  field  expansions  are  substituted  into 
Maxwell's  equations,  yielding  [Ref.  7],  after  judicious  cross 
substitutions  the  following  modal  field  generating  equations  in  terms 


of  two  scalar  potentials,  R,Z,m  )  and  r2(  R,Z,m  ): 

<J>  xem  =  j  4  (  m<(>xVIj  -  Rpr  Vr2  )  (A.5) 

♦  •  em  =  S-  (A.6) 

♦  xhm  =  j  fm  (  m^xVP,  +  Rervr,  )  (A.7) 

♦  ■  h„  -  £  (A.8) 

where  the  two  dimensional  gradient  operator  is  defined  as 

V  =  [  R  <d/dr >  +  Z  (a/az)  ]  (A. 9) 

and  fm  is  a  multiplicative  function  given  by 

UR.Z)  =  [  Hr(R.2)  ^ (R.Z5  R2- m2  j-1  .  (A.10) 

Note  that  the  potentials,  1^  and  r2  are  proportional  to  the  azimuthal 


components  of  the  electric  and  magnetic  fields,  em  and  hm.  These 

potential  functions  are  continuous  everywhere,  including  at  dielectric 
and  magnetic  interfaces.  This  property  of  uniform  continuity  in  the 
field  is  used  to  obtain  a  numerical  solution  of  the  fields. 
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The  total  fields  may  be  described  in  term  of  the  uniformly 
continuous  scalar  potentials  rj(  R,Z,m  )  and  r2(  R,Z,m  )  such  that 

em(R,Z)  =  j  fm  (mVIl  -  Rm*xvr2)  +  *2-  (A.11) 

and 

hm(R,Z)  =  j  f„  (mvr2  +  Rer»xVi;)  +  ♦  ^  .  (A.l  2) 

This  method  of  expanding  the  electromagnetic  fields  into  modal 
components  is  named  the  Coupled  Azimuthal  Potential  (CAP)  method 
due  to  the  fact  that  the  scalar  potentials  T^R,  Z,  m)  and  r2(R,  Z,  m) 
are  coupled  in  these  field  equations  and  are  not  separable. 

Due  to  the  symmetry  in  the  azimuthal  coordinate,  <J>,  in  radome 
design,  an  arbitrary  meridional  cross  section  of  the  radome 
construction  may  be  represented  in  an  (R,  Z)-pIane  such  that  there 
exists  a  two-dimensional  solution  domain  S,  with  a  boundary  curve 
forming  the  surface  of  revolution  5S.  It  has  been  shown  by  Morgan, 
Chang  and  Mei  [Ref.  7]  that  in  this  solution  domain  the  potentials,  T j 
and  T 2,  satisfy  the  following  formally  self-adjoint  pair  of  partial 
differential  equations  (PDE's), 

V(fm  Rervn)  +  -jfr,  +  mv-(k*xvr2)  =  0  (A.l  3) 

and 

-mv-i^xvn)  +  V ■  C ^  rm,  vr2)+-|-r2  =  o  ,  (A.i4) 

with  Dirichlet  boundary  conditions  (BC)  for  Tj  and  r2  specified  on  5 S  . 
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In  addition  to  this  PDE  formulation  there  is  an  easily  managed 
variational  principle  [Ref.  7]  established  using  either  an 
Euler-Lagrange  approach  or  a  stationary-theorem  approach.  The 
variational  principle  is  based  upon  the  existence  of  a  stationary 
functional  of  the  potentials,  1^  and  r2,  and  their  gradients  such  that 

sf  =  jj  L(R,zr1,r2,vr1,vr2)  dRdz  ,  (a.15) 

s 

where  the  Lagrangian  operator,  L,  is  unique  within  an  arbitrary 
constant  multiplier  and  arbitrary  independent  additive  function.  The 
Lagrangian  operator  has  the  quadratic  form 

L  =  fjn  VTj-  (RejVrj  +  m4>xVr2) 

+  vr2‘  (  Rur  vr2  -  m0X  VT]  ) 

-  i  (er  if  +  M,  r22  )  .  (A.16) 

With  this  analytical  approach  to  the  solution  of  the  electromagnetic 
field's  behavior  through  the  axially  symmetric  radome  structure,  the 
determination  of  the  field  on  a  point-to-point  basis  is  accomplished 
through  a  numerical  finite-element  algorithm. 

B.  FINITE  ELEMENT  ALGORITHM 

The  numerical  solution  of  the  electromagnetic  field  problem 

involving  scattering  by,  and  penetration  of,  the  radome  is 

accomplished  by  the  development  of  a  finite-element  algorithm 

using  the  mesh  geometry  in  Fig  A. 2.  This  configuration  in  the 

meridional  (R,  Z)-plane  uses  linear  triangular  elements  which  are 
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Unbounded 
Exterior 
REGION  III 


Figure  A. 2  Finite  Element  Mesh  Geometry 
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arranged  by  the  EMRAD  program  to  conform  to  the  surfaces  of  the 
axisymmetric  radome.  The  selection  of  surface  conformation 
simplifies  numerical  analysis  of  the  normal  derivative  discontinuities 
at  the  surface  boundaries.  The  two-dimensional  mesh  configuration 
that  is  developed  is  named  the  semiannular  conformal  (SAC)  mesh. 

Figure  A. 2  graphically  represents  the  separate  spatial  regions 
which  are  of  interest.  The  homogeneous  inner  spherical  core  of  the 
radome  occupies  Region  I.  The  inhomogeneous  axisymmetric  radome 
wall  construction  cross-section  is  in  Region  II.  The  unbounded 
exterior  region  exists  exterior  to  the  outer  sphere  in  the 
homogeneous  free-space  of  Region  III.  Note  that  the  surfaces  of  the 
inner  sphere,  Region  I,  and  the  exterior  Region  III,  bound  the  area  of 
interest  for  numerical  analysis  in  Region  II. 

The  electromagnetic  fields  within  the  homogeneous  inner  core 
of  the  radome  (Region  I)  and  the  unbounded  exterior  (Region  III)  are 
defined  in  terms  of  the  appropriate  spherical  harmonic  modal 
expansions  of  the  radial  TE  and  TM  field  components.  At  the 
surfaces,  r  =  a  and  r  =  b,  the  modal  elements  in  these  fields  are 
applied  as  sets  of  Dirichlet  boundary  conditions.  Then,  the  numerical 
solution  for  the  fields  in  the  finite-element  domain  of  Region  II  is 
generated  for  each  boundary  condition.  Finally,  the  unknown  modal 
expansion  coefficients  for  the  core  and  exterior  scattered  fields  are 
obtained  by  enforcing  continuity  of  the  electromagnetic  fields  at  the 
interfaces  along  r  =  rj  and  r  =  r2.  This  approach  to  numerical  solution 
of  the  fields  is  the  basis  of  the  unimoment  method.  In  the 
implementation  of  the  variational  principle  using  the  finite-element 


method,  the  azimuthal  modal  field  components  are  selected  as  the 
fundamental  unknowns. 

In  Region  II  of  Fig.  A. 2,  the  SAC  mesh  is  composed  of  numerous 
interior  nodes.  Each  of  these  interior  nodes  is  surrounded  by  six 
triangular  regions  known  as  "elements".  The  azimuthal  field 
components  within  the  mesh  are  represented  as  basis  function 
expansions, 

N 

$  •  em(R,Z)  =  £  em(n)  un(R,Z)  (A. 17) 

n=l 

and 

♦  •  Hm(R,Z)  =  f  hm(n)  u„(R,Z)  .  (A.l  8) 

n=l 

Each  of  the  respective  sets  of  em(n)  and  hm(n)  represent  one  of  the  N 
complex  nodal  values  found  in  the  azimuthal  fields.  The  finite 
element  basis  functions  un(R,  Z)  are  area  coordinate  linear  pyramid 

functions  which  each  assume  a  value  of  unity  at  their  associated 
node,  n,  and  zero  at  each  of  the  neighboring  nodes.  There  is  an 
assumption  of  linear  variation  within  each  of  the  triangular  elements 
in  the  surrounding  group  such  that 

un(R,Z)  =  an  fR  +  bn j  Z  +  c^f  .  (A.l  9) 

The  real  constants  (a,  b,  c)  corresponding  to  each  node,  n,  in  element 

C  is  found  in  terms  of  the  nodal  (R,  Z)-coordinates  of  the  element. 

The  variational  finite-element  technique  is  initiated  by 

substituting  the  linear  shape-function  expansions  in  Equations  (A.  17) 

and  (A. 18)  into  the  stationary  Euier-Lagrange  functional.  SF.  as 
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described  in  Equation  (A.  15).  The  stationary  point  of  SF  is  then 
computed  by  imposing  the  null  derivative  condition  with  respect  to 
each  of  the  unknown  coordinates  representing  the  interior  nodal 
fields,  while  constraining  the  "known"  boundary  nodal  fields.  The 
resulting  linear  system  describing  the  nodal  fields  is 

N  r  r 

I  (em(i)JJ  Um  RerV(Ru„)  V(RUi)-  Refu„Ui]dRdZ 

i=1  Sn 

+hm(i)JJ  [mfm  V(Run)-  4>x  V( RUj) ]dRdZ )-0  (A. 20) 

Sn 

and 

£(MOjJ  tt„  RnrV(RUn)  V(Rui)-R(M.nUjldRdZ 
sn 

[m^V(Run) .*x  V(RUi)]dRdZ)=0  (A.2  1) 

s„ 

for  each  value  of  the  internal  nodes,  n,  where  i  spans  all  nodes, 
including  boundary  nodes.  As  these  above  equations  are  enforced 
for  each  internal  node,  they  involve  integrations  over  only  the 
elements  connected  to  this  node,  Sn,  and  therefore  relate  the  value  of 
em(n)  and  hm(n)  oniy  to  adjacent  nodal  values,  em(i)  and  hm(i).  Then, 

in  the  SAC  mesh  implementation,  in  the  case  of  each  internal  node 
and  the  six  adjacent  nodes,  the  above  equation  will  relate  14 
adjacent  node  azimuthal  field  values,  some  of  which  may  be  specified 
at  the  boundary. 
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The  computation  for  each  nodal  value,  as  indicated  by 
Equations  (A.20)  and  (A.21),  results  in  a  sparse  global  system  matrix 
which  contains  generally  complex  elements.  At  each  node  the  mean 
values  of  the  spatially  variable  er(R,  Z)  and  pr(R,  Z)  are  used  as  the 
material  constants  in  that  triangular  element.  Applying  the  Dirichlet 
boundary  conditions  to  the  SAC  mesh  shown  in  Fig.  A. 2,  Morgan  and 
Mei  [Ref.  5]  found  the  numerical  finite-element  solution  for  the 
electromagnetic  fields  in  Region  II  using  a  two-sweep  block-by-block 
matrix  inversion  algorithm  related  to  the  Riccati  transform  for 
difference  equations.  At  each  constant  colatitude,  0  =  0k,  a  vector  of 
internal  nodal  azimuthal  field  values  along  the  ray  from  r}  to  r2  is 
defined  by  P  k  .  The  resultant  global  system  matrix  defined  by 
Equations  (A.20)  and  (A.21)  has  a  diagonal  tri-block  submatrix 
structure,  which  can  be  represented  locally  by  relationships  between 
adjacent  radial  node  vectors  in  the  mesh, 

Ak  Pk-1  +  Pk  +  Q  Pk+1  =  ak  *  (A. 2  2) 

such  that  Ak,  Bk,  and  Ck  are  sparse  submatrices.  These  form  the 
tri-block  diagonal  structure  of  the  global  system  matrix  while  a  k  , 
the  local  driving  function,  is  generated  from  Equations  (A.20)  and 
(A.21)  using  the  known  Dirichlet  BC's  at  r=a  and  r=b. 

The  system  solution  in  the  two-dimensional  meridional  (R,  Z) 

plane  is  performed  in  a  step-wise  manner  beginning  at  the  positive 

Z-axis  (0k  =0]  =0)  and  proceeding  in  equal  marching  increments  of 

A0.  At  each  0k  ,  the  SAC  mesh  is  constructed  using  a  minimum 

distortion  scheme  to  conform  to  the  approximate  piecewise  smooth 
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surfaces.  At  each  0k  increment,  as  the  mesh  is  constructed  in  this 
sequential  circumferential  manner,  the  three  block  matrices  are 
generated  which  relate  adjacent  radial  node  vectors. 

During  the  incremental  marching  of  0k  ,  the  block-by-block 
Riccati  matrix  inversion  is  implemented.  The  Riccati  matrix,  Rk  ,  and 
the  associated  vector,  Sk  ,  are  related  to  the  vector,  P  k  ,  by  the 
transformation 

ftc-l  =  Rk  Pk  +  Sk  .  (A. 2  3) 

Substituting  this  Riccati  relation  into  Equation  (A. 22),  solving  for  P  k 
and  returning  to  the  form  of  the  Riccati  transform  yields  the 
recurrence  relations 

Rk+i  =  — (  Bfc  +  Ak»Rk  )  '•Q  (A. 24) 

and 

Sk+1  =  (  ®k  +  ^k«Rk  )  1#(  Pk  “  Ak»Sk  )  •  (A. 2  5) 

Using  these  relations,  the  Riccati  transform  sweeps  forward 
from  initial  conditions,  using  the  equations  of  step-wise  evolution  in 
Equations  (A. 24)  and  (A. 25)  to  sequentially  generate  and  store  the  Rk 
matrices  and  the  Sk  vectors.  This  is  called  the  "first  sweep”.  The 
solution  vectors  P  k  are  then  obtained  by  using  the  transform  in 
Equation  (A. 23)  to  perform  a  "backsweep"  with  calculations  using  the 
stored  Riccati  arrays.  This  sequentially  localized  technique  for 
handling  the  Dirichlet  boundary  value  problem  blends  easily  with 
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the  unimoment  method  which  is  used  to  solve  the  electromagnetic 
field  in  Regions  I  and  III  of  Fig.  A. 2. 

C.  UNIMOMENT  METHOD 

The  unimoment  method  is  a  computational  procedure  which 
provides  for  the  numerical  solution  of  the  incident,  scattered  and 
spherical  core  electromagnetic  field  expansions  when  dealing  with 
the  constructed  SAC  finite-element  mesh  of  Fig.  A. 2  The  numerical 
solution  is  generated  by  representing  the  core  (r<rj)  and  the  exterior 

scattered  (r>r2)  fields  in  terms  of  radial  TE  and  TM  spherical 
harmonic  expansions,  each  with  unknown  core  and  scattered  modal 
field  coefficients.  In  order  to  solve  for  these  unknown  coefficients, 
the  core  expansion  of  Region  I  is  applied  as  a  boundary  condition 
along  r  =  a,  while  the  boundary  condition  at  r  =  b  is  determined  as 
the  addition  of  the  incident  field  and  the  scattered  field  expansion  in 
Region  III. 

Superposition  is  then  applied  to  the  boundary  conditions, 

composed  of  the  core  field  expansion  and  the  exterior  field  expansion 

in  order  to  find  a  complete  solution.  The  complete  field  solution  is 

found  from  the  addition  of  the  partial  field  solutions  generated  from 

two  sets  of  applied  BC's.  The  first  set  of  boundary  conditions  applied 

is  composed  of  the  modal  spherical  harmonic  core  fields  along  r  =  a 

with  zero  fields  along  r  =  b.  The  second  set  of  boundary  conditions 

contains  the  spherical  harmonic  exterior  scattered  and  incident  fields 

at  r  =  b  with  null  fields  at  r  =  a.  These  multiple  boundary  conditions 

generate  multiple  solutions  at  the  nodes  along  both  boundaries, 
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r  =  rj  and  r  =  r2.  The  unknown  core  and  scattered  modal  field 

coefficients  are  obtained  by  setting  the  weighted  superposition  of 
interior  solutions  equal  to  the  original  field  expansions.  This  equality 
is  enforced  along  the  boundaries  at  f  =  ri  and  r  =  r2  ,  utilizing  least 
squares. 

The  success  of  the  unimoment  method  is  contingent  upon  the 
capability  to  analytically  represent  the  electromagnetic  fields 
everywhere  outside  of  the  finite  element  mesh  region.  These 
representations,  as  developed  for  Regions  I  and  III,  are  given  as 
infinite  series  of  spherical  harmonics  which  must  be  truncated  to 
facilitate  numerical  calculations.  The  plane-wave  fields  that  are 
incident  on  the  radome  are  decomposed  in  spherical  coordinates  into 
azimuthal  modes  as 


M 


E?A(r,fy|>)  S  £  em,0(r’e) 


m=0 


and 


M 


ri0tf6(r,a<t>)  =  £  hm,o(r’0) 


m=0 


jsin(m<{>) 

[_  cos(mb)J 

cos  (mb) 

L  jsin(mb)J 


(A. 26) 


(A. 27) 


where  M  is  the  mathematical  truncation  index  of  the  series 
representation.  The  bracketed  sinusoids  represent  the  polarization 
convention  used  for  TE  and  TM  modes,  respectively. 

During  computations,  the  numerical  solution  of  the  unknown 


core  and  scattered  modal  field  coefficients  is  run  concurrently  with 

the  sequential  Riccati  transform  finite-element  solution  in  Region  II. 
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On  the  forward  sweep,  boundary  conditions  for  both  the  incident  and 
the  modal  fields  are  generated  only  at  the  local  node  row  being 
considered  in  that  particular  sweep  step.  On  the  backsweep  the 
nodal  values  for  each  applied  boundary  condition  are  obtained  along 
the  r  =  Tj  and  r  =  r2  contours.  The  corresponding  incident  field  and 

spherical  harmonic  field  at  each  node  are  generated,  and  then  the 
nodal  residual  fields  are  formed.  During  the  backsweep,  inner 
products  integrations  are  accumulated  and  their  contributions  are 
added  to  a  field  moment  matrix  and  the  incident-field  driving  vector, 
which  corresponds  to  each  incident  field.  The  field  expansion 
coefficients  for  each  incident  field  are  obtained  by  multiplying  the 
inverse  of  the  field  moment  matrix  with  that  field’s  driving  vector. 
With  the  field  expansion  coefficients  calculated  and  collected,  the 
electromagnetic  field  within  the  radome  coTe  can  be  determined  on  a 
point-to-point  basis  as  developed  in  Chapter  II  . 
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Appendix  B 


RADSHP  PROGRAM  CODE 


PROGRAM  RADSHP 
C 

Q  ************************************************************** 

c 

C  GENERATING  THE  MATERIAL  FILE  AND  DATA  STRUCTURE  FILE 
C  FOR  A  LAYERED  DIELECTRIC  RADOME 
C  MATERIAL  FILE  PREPARED  FOR  EMRAD 
C  DATAFILE  PREPARED  FOR  CURVE  DIGITIZER 
C  USE  EMCADIN  TO  GENERATE  STRUCTURE  FILE  FOR  EMRAD 
C  Written  by  R.J.Vince  Nov-Dec  89 
C 

£  ************************************************************** 

C  VARIABLE  DECLARATIONS 

Q  *************************************************************** 

REAL  X(5,1001),  Z(5,1001),  PI,  DELTHETA,  DEGTORAD 
REAL  THETA,C(5),R(5),H(5),DELZ(5),SHFTFACTORF(5,7) 

INTEGER  IBIG,  I,  L,  SEPRATN,  NSHAPE,  NLMAX,  RESOLV,  RESPTS 
COMPLEX  ER(5),UR(5) 

CHARACTER'  10  SHAPE(6) 

CHARACTER‘64  HDRI.HDR2 
CHARACTER  *  8  YNAME 
CHARACTER*  12  MATFIL.DATFIL 
CHARACTER*!  YN 


SLARGE 

C 

Q  *************************************************************** 

C  INITIAL  VALUES 

p  *************************************************************** 


no 


PI=3. 14 15927 
DEGTORAD  =  Pl/180. 
SEPRATN  =  999990 
RESOL  V=  1000 
RESPTS=RESOLV+ 1 
NLMAX  =  5 


C  mTIALIZE  SHAPE  FUNCTION  FACTOR  ARRAY 
C 

DO  67  1=1  NIL  MAX 
DO 65  J=l,6 
F(IJ)=0.0 
65  CONTINUE 
67  CONTINUE 
C 

C 

£  ****************************************************************** 

C  INPUT  NECESSARY  DATA 

£  ****************************************************************** 

C 

C  PROMPT  FOR  OUTPUT  FILENAMES 
C 

WRITER,  110) 

1 10  FORMAT(//9X, 'WELCOME  TO  RADSHP’, 

1  ///9X.THIS  PROGRAM  GENERATES  RADOME  MATERIAL  AND  SHAPE  DATA', 

2  //7X/FILES  FOR  USE  WITH  EMRAD,  BASED  ON  USER  INPUT.', 

3  ////9X.TWO  FILES  ARE  REQUIRED  TO  HOLD  PROGRAM  OUPUT  DATA.', 

4  I//9X,'  PLEASE  ENTER  THE  NAME  OF  THE  RLE  TO  BE  USED  TO  ', 

5  /7X/HOLD  THE  PROGRAM  MATERIAL  PARAMETER  DATA.  THE  EXTENSION', 

6  /7X,'.MAT  WILL  BE  ADDED  AUTOMATICAI  LY,  I.  E.  YOURNAME.MAT',/) 
READ(*,100)  YNAME 

NCHAR  =  INDEX(YNAME,' ')  -1 
MATFIL  =  YNAME(LNCHAR)  //  ’.MAT 

OPEN(  1  pILE=MATFIL,STATUS='UNKNOWN',ACCESS=’SEQUENTIAL'. 

1  1  1 


1  FORM=TORMATTED') 

WRITE(Mll) 

1 1 1  FORMAT(///9X,’PLEASE  ENTER  THE  NAME  OF  THE  FILE  TO  BE  USED 

1  /7X.TO  HOLD  THE  PROGRAM  OUTPUT  DATA.  THE  EXTENSION  .DAT  WILL', 

2  /TX.’WELL  BE  ADDED  AUTOMATICALLY,  I.  E.  YOURNAME.DAT  J1X) 
READ<*,100)  YNAME 

NCHAR  =  INDEX(YNAME,’  ■)  -1 
DATFIL  =  YNAME(1  :NCHAR)//  '.DAT 

OPEN(2JILE=DATFIL,STATUS=,UNKNOWN’^CCESS=,SE0UENTIAL,> 

1  FORM=EORMATTED) 


C 

C  PROMPT  FOR  #  OF  LAYERS  WITHIN  THE  RADOME 
C 

5  WRITE(\114)NLMAX 

1 14  FORMAT(//7X,ENTER  THE  NUMBER  OF  LAYERS  (.LE.  ,I2/)7/7X) 
READ/*,*)  NL 

IF  ( (NL  JLT.  1)  .OR.  (NL  ,GT.  NLMAX) )  THEN 
WRITE(* ,  1 1 5)NLMAX+ 1 

115  FORMAT(//7X,'NUMBER  OF  LAYERS  IS  OUT  OF  RANGE,(Ck#<M2,')7/7X) 
GOTO  5 

END  IF 


C  LOOP  TO  INPUT  LAYER  INFORMATION 
C 

WRITE(*,7) 

7  FORMAT(/8X,Each  layer  of  the  radome  will  have  material  consts,', 

1  / lOX.'an  outer  height  of  the  radome  layer.H,’, 

2  /10X,'a  outer  base  radius.R,  and', 

3  /10X,'a  shape  function,  &  shape  factor  parameters  (if  appl.),', 

4  //I OX, The  first  "layer"  will  be  the  inner  core  of  the  radome,  ’, 

5  /10X,'normally  filled  with  air.’) 
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DO  1 2  L=  1  ,NL 
39  CONTINUE 

WRITE(*,9)  L 

9  FORMAT(/,'l  Layer  No:  '  J3) 

IF  (L.EQ.l)  THEN 

WRITE(*,*)  ’  This  is  the  inner  core  of  the  radome' 
ELSEIF  (L.EQ.NL)  THEN 

WRITE(*,*) '  This  is  the  outer  layer  of  the  radome’ 
END  IF 


C  Prompt  for  Height  and  Radius  of  this  layer 
C 

WRITE(V) 

WRITER  ,43) 

43  FORMAT(/8X/Please  enter  the  radome  layer  height,H,’ 
1  /8X,'  and  the  outer  base  radiusJR,  for  this  layer.’) 

READf*,*)  H(L),R(L) 


C  Prompt  for  shaping  function,  to  determine  equations  to  be  used 
C  to  generate  this  layer  of  the  radome 
C 


WRITE(V) 

WRITE(*,10) 

10  FORMAT(/8X,The  following  shaping  functions  are  available,  ’ 

1  /8X,’  1.  General  Ogive' 

2  /8X,'  2.  Tangent  Ogive' 

3  /8X,'  3.  von  Karman’ 

4  /8X,'  4.  Power  Series' 

5  /8X,'  5.  Parabola,  curved' 

6  /8X,'  6.  Parabola,  capped  by  pointed  apex' 

7  //8X,'  Indicate  desired  shaping  function  by  entering  a  number' 

8  /10X,’  1,2, ...  or  6  ’) 
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READ<*  ,*)  NSHAPE 


IF  (NSHAPE.EQ.1)  THEN 
Ogive 

SHAPE(L)  =  'Gen.  Ogive' 

CALLOGIVE(L,X,Z,H,R,F.RESOT  V.RESPTS.NLMAX) 

ELSEIF  (NSHAPE.EQ.2)  THEN 
Tangent  Ogive 

SHAPE(L)  =  TAN  Ogive  ’ 

CALL  TANOGIVE(L,X,Z,H,R  P.RESOLV  .RESPTS.NLMAX) 

ELSEIF  (NSHAPE.EQ.3)  THEN 
von  Kami  an 

SHAPE(L)  =  von  Karman’ 

CALL  VONKRMN(L,X2,H^JlESOLV>RESPTS>NLMAX) 

ELSEIF  (NSHAPE.EQ.4)  THEN 
Power  Series 

SHAPE(L)  =  Pwr  Series' 

CALL  POWER(L,X,Z,H,RTF,RESOLV,RESPTS  .NLMAX) 

ELSEIF  (NSHAPE.EQ.5)  THEN 
Parabola  (curved) 

SHAPE(L)  =  Para  Curve' 

CALL  PARACRV(L,X,Z,H,R,F,RESOLV,RESPTS,NLMAX) 

ELSEIF  (NSHAPE.EQ.6)  THEN 
Parabola  with  apex 

SHAPE(L)  =  Para  APEX' 

CALL  APEX(LX,Z,H,R>F,RESOLV,RESPTS  .NLMAX) 

ELSE 

Selection  out  of  limits 

WRITE  (*,*)  ’Selection  is  not  available,  try  again.’ 

WRITE  (* ,*) 
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GOTO  39 
END  IF 


C  Provide  listing  of  inner  layer  parameters 
C 

WRITE(*,554) 

DO  11  1=1  J- 

WRITE  (*,555)  I,SHAPE(I),H(I)^(I),(F(U)J=1,7) 

11  continue 

WRITE(*  *) 

12  CONTINUE 

C 

C  SHIFT  RADOME  SO  THAT  ORIGIN  IS  NEAR  THE  CENTER  OF  INNER  LAYER 
C  NOTE:  FACTOR  is  adjusted  for  optimal  placement  of  origin  for 
C  mesh  construction  by  EMRAD 
C 

FACTOR=0.4 


C  Max  inner  height  =  Z(NLJRESPTS) 

C 

SHFT=Z(1  ,RESPTS)*(F  ACTOR) 

DO  88  L=  1  ,NL 
DO  77  I=I,RESPTS 
Z(L,I)=Z(LJ)-SHFT 
77  CONTINUE 

88  CONTINUE 


C 

Q  ****************************************************************** 

C  OUTPUT 

£  ****************************************************************** 

115 


DO330L=ltNL 

C  OUTPUT  LAYER  XX  POINTS  FROM  TIP  TO  X-Z  PLANE 
DO  220  I=RESPTS,1,-1 
WRITE(2 ,*)  X(L,I),Z(L,I) 

220  CONTINUE 

C  CLOSE  LAYER  BELOW  SHIFTED  X-Y  PLANE 
C  AT  Z=[  -1*[SUM  OF  (LAYER  THICKNESS)]  -  SHFT  ] 

C 

c 

C  DETERMINE  LAYER  THICKNESS,  INNER  LAYER  HAS  NO  SHIFT  FOR  THICKNESS 
C  OTHER  LAYERS  HAVE  BOTTOM  THICKNESS  DETERMINED  BY  RADIAL  THICKNESS 
C  NOTE  -  THE  SHIFT  PERFORMED  FOR  THICKNESS  IS  ADDITIVE  FOR  LAYERS 
C 

IF  (L.EQ.l)  THEN 
THICK  =  0 
ELSE 

THICK  =  THICK  +  R(L)-R(L-1) 

ENDIF 

C 

BOTTM  =  (-1.0  *  THICK)  -  SHFT 
WRITE(2 ,*)  R(L),  BOTTM 
WRITE(2,*)  0.0,  BOTTM 
WRITE(2,*)  SEPRATN.SEPRATN 
330  CONTINUE 

C  SIGNAL  TERMINATION  OF  DATA  SET 
WRITE(2,100)  END,  END' 

C  Put  summary  of  data  to  screen 
C 

WRITE(*,554) 

DO  440  1=1  NL 

WRITE  (*,555)  I,SHAPE(I),H(I).R(I),(F(U)J=1,7) 
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440  CONTINUE 


C 

C  PROMPT  FOR  HEADER  INFORMATION 
C 

WRITE(M12) 

1 12  FORMAT(//7X,THE  HEADERS  ALLOW  THE  USER  TO  IDENTIFY  THIS7/7X, 
l'SET  OF  DATA  FROM  ALL  OTHER  SETS.’, 

2//7X, 'PLEASE  ENTER  HEADER  #1  (64  CHARACTERS  MAX)' J/1X) 
READ(*,100)  HDR1 
WRITER, 113) 

1 13  FORMAT(//7X/PLEASE  ENTER  HEADER  #2  (64  CHARACTERS  MAX)7/7X) 
READ(MOO)  HDR2 


C 

C  WRITE  HEADER  INFORMATION  INTO  MATFRIAL  OUTPUT  DATA  FILE 
C 

WRITE(  1,105)  HDR1 
WRITE(  1,105)  HDR2 

C 

C  WRITE  LAYER  INFORMATION  INTO  MATERIAL  OUTPUT  DATA  FILE 
C 

WRITE(1 ,106)  NL 

C  Prompt  for  material  Er  for  each  layer  of  construction 
C 

D0  812L=1JML 
WRITE(*,9)  L 

IF  (L.EQ.l)  THEN 

WRITE(*,*) '  This  is  the  inner  core  of  the  radome' 

ELSEIF  (L.EQ.NL)  THEN 

WRITE(*,*) '  This  is  the  outer  layer  of  the  radome' 

ENDIF 
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C  Prompt  for  material  constants 
C 

WRITE(*,*)  'Enter  Er  for  this  Layer: 1 
READCV)  ERR 

C  Ensure  that  Imag  Er  is  <  0 
C 

ERI  = -0.000001 

ER(L)=CMPLX(ERR£RI) 

UR  (L)=C  MPLX  ( 1 .0,- 1  .E-6) 
WRITE(1,102)  ER(L),UR(L) 

812  CONTINUE 


C  PLACE  COPY  OF  CONSTRUCTION  AT  BOTTOM  OF  MATERIAL  FILE 
C 

WRITE(  1,554) 

DO  660  1=1  M- 

WRITE  (1.555)  I,SHAPE(I),H(I)^(I).(F(U)J=1,7) 

660  CONTINUE 


C  —  FINAL  OUTPUT  TO  SCREEN 

WRITE(*,999)  MATFIL.DATFIL.DATFIL 
999  FORMAT(//8X, 

L  Computations  and  output  are  completed, the  datafile  ,  ,A, 

2/.8X,’  holds  the  material  specifications  for  input  to  EMCAD.’ 
3//.8X,’  The  datafile,  ’.A,’,  holds  the  radome  shape  data  for’ 

4/,8X,’  viewing  in  the  CAD  package,  CURVE  DIGITIZER.  ’ 
5//.8X,'  The  structure  file  required  for  input  to  EMCAD  can  be' 
6/,8X,'  generated  by  the  executing  the  program  EMCADIN,' 
7/.8X,'  using,  ',A,' ,  as  input  to  EMCADIN, ' 

8//,8X,'  The  output  file  from  EMCADIN  will  be  in  proper  format' 
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9/.8X,'  for  input  to  EMCAD  for  the  structure  file.') 

WRITE(V) 

C 

Q  ****************************************************************** 

100  FORMAT(A) 

101  FORMAT(I5) 

102  FORMAT(4(E14.6)) 

103  FORMAT(I5,5(E14.6)) 

104  FORMAT(2(I5)) 

105  FORMATf  » 

106  FORMAT(I5) 

554  FORMAT  (///////////////// 

1  //lX/LAYER'^X,'  SHAPE  ',2X, 'HEIGHT', 2X, 'RADIUS', 

1  2X,'  FI  ’,2X,'  F2  ',2X,'  F3  ',2X,'  F4  ’,2X,'  F5  ’, 

2  2X,’  F6  ',2X,'  F7  71 X.' . ',2X,' - ’,2X,‘ - ', 

3  2X,' . ',2X, 2X,'-~’,2X,’—',2X,'-~’, 

4  2X,’ — ',2X,' — ’,2X,’ — ') 

555  FORMAT  (1X,I5,2X,A10,2X,F6.2,2X,F6.2,7(F6.2)) 

9999  CONTINUE 
C 

C  END  PROGRAM  RADOME 
C 

CLOSE(l) 

CLOSE(2) 

STOP 

END 

********** ************ ************* 

C  OGIVE  SHAPING  FUNCTION  SUBROUTINE 

SUBROUTINE  OGIVE(L,X,Z,H  ,R  F,RESOLV,RESPTS,NLMAX) 

INTEGER  IBIG,  I,  L,  SEPRATN,  NSHAPE,  NLMAX,  RESPTS,  RESOLV 
REAL  X(NLMAX .RESPTS),  Z(NLMAX  JIESPTS),  PI,  DELTHETA,  DEGTORAD 
REAL  THETA,C(5)JR(5),H(5),DELZ(5),SHFT .FACTOR 
REAL  ALPHA  1(5),  ALPHA(5),  B,  M,  RL, GAMMA,  F(NLMAX,7) 

COMPLEX  ER,UR 

1  1  9 


CHARACTER*!  YN 


C  SET  CONSTANTS 
C 

PI=3. 14 15927 
DEGTORAD  =  PI/180. 

C 

C  DEFINE  START  AND  END  POINT  OF  EACH  RADOME  LAYER 
C 

Z(L,1)=0.0 

X(L,1)=R(L) 

Z(L,RESPTS)=H(L) 

X/L,RESFTS)=0.0 

C . DEFINE  RADOME  LAYER  BY  OGIVE  EQUATIONS  .... 

DELZ(L)=H(L)/RESOLV 

C  Prompt  for  Shaping  factors  for  this  ogive  layer 
C 

WRITE/*,*) 

WRITE(*,43) 

43  FORMAT(8X,'The  ogive  shaping  function  requires  three  shaping’, 
l/8X,'factors  for  determining  the  radome  construction.' 
l//8X,Thc  first  factor,  FI ,  determine  the  arc  of  curvature', 
l/8X,’of  the  radome  surface.  Please  enter  FI  :’) 

REAL>(*  *)F(L,1) 

WRITE/*,*) 

WRITE/*  ,45) 

45  FORMAT/8X,'  Two  other  shaping  factors  are  used  to  determine  the', 
l/8X,’relative  position  of  the  center  of  the  ogive  shaping', 
l/8X,'arc,  relative  to  the  radome  apex.’, 

1//8X,'  The  first  of  these  factors,  F2,  determines  the  perp.  ', 
l/8X,'distance  from  the  z-axis  to  the  center  of  the  arc.', 
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1//8X,'  The  second  positioning  factor,  F3,  determines  the', 
l/8X,'the  distance  from  the  center  of  the  arc  to  the  apex', 
l/8X,'of  the  radome,  along  the  z-axis.', 

1//8X, Please  enter  the  position  factors  for  the  center  of  the’, 

1/8X,'  curvature  arc,  F2  and  F3.1) 

READ(*,*)F(L,2),F(L,3) 

DO  33  I=2,RESPTS-1 

Z(L,I)=Z(L,1- 1  )+DELZ(L) 

C  EQUATION  OF  RADOME  IS  x=SQRT(  Fl**2-(F2-(H-z))**2 )  -  F3 
C 

X(LJ)  =  SQRT(  F(L,1)**2  -( F(L,2)  -  (H(L)-Z(LJ))  )**2  )-F(L,3) 

33  CONTINUE 


C . END  OF  RADOME  LAYER  .... 

RETURN 

END 

t*********************!^^  Qp  QQjyp 
*********************************** 

C  TANGENT  OGIVE  SHAPING  FUNCTION  SUBROUTINE 

SUBROUTINE  TANOGrVE(L,X,Z,H,R,F,RESOLV, RES  PTS.NLMAX) 

INTEGER  IBIG,  I,  L,  SEPRATN,  NSHAPE,  NLMAX,  RESPTS,  RESOLV 
REAL  X(NLMAX .RESPTS),  Z(NLMAX PESPTS),  PI,  DELTHETA,  DEGTORAD 
REAL  THETA, C(5)P(5),H(5)J)ELZ(5),SHFTPACTOR 
REAL  ALPHA  1(5),  ALPHA(5),  B,  M,  RL, GAMMA,  F(NLMAX,7) 

COMPLEX  ER,UR 
CHARACTER*  1  YN 

C  SET  CONSTANTS 
C 

PI=3. 1415927 
DEGTORAD  =  PI/ 180. 
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c 

C  DEFINE  START  AND  END  POINT  OF  EACH  RADOME  LAYER 
C 

Z(L,1)=0.0 

X(L,1)=R(L) 

Z(L,RESPTS)=H(L) 

X(L,RESPTS)=0.0 

C . DEFINE  RADOME  LAYER  BY  TANGENT  OGIVE  EQUATIONS  .... 

DELZ(L)=H(L)/RESOLV 

C  Determine  Shaping  factors  for  this  ogive  layer 
C 

F(L,4)  =  R(L)/2.0  +  2.0*(  H(L)**2,0  )/R(L) 

F(L,5)  =  F(L,4)  -  R(L)/4.0 

WRTTE(*  ,*) 

WRITE(*,43)F(L,4),F(L,5) 

43  FORMAT (8X .The  tangent  ogive  shaping  function  determines  two', 
l/8X,'factors  for  the  radome  construction.  These  are  calculated' 

1/8X, "based  on  the  height  and  radius  of  the  radome  layer.', 

1/8X, These  factors  are  F4  =  ',F6.2,'and  F5  =  ’,F6.2) 


D0  33I=2,RESPTS-1 

Z(L,I)=Z(L  ,1- 1  )+DELZ(L) 

C  EQUATION  OF  RADOME  IS  x=SQRT(  F4**2  -  z**2 )  -  F5 
C 

X(LJ)  =  4.0*  (SQRT(  F(L,4)**2-  Z(LJ)**2  )  -  F(L,5)) 

33  CONTINUE 


C . END  OF  RADOME  LAYER  .... 
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RETURN 

END 

*»»**»»***»******»**»*END  op  JAN  OGIVE 
••*«*•»*********»•*• •*••••••••*»*•• 

C  von  Karman  SHAPING  FUNCTION  SUBROUTINE 

SUBROUTINE  VONKRMN(L,X,Z,H,R,RESOLV,RESPTS,NLMAX) 

INTEGER  IBIG,  I,  L,  SEPRATN,  NSHAPE,  NLMAX,  RESPTS,  RESOLV 

REAL  X (NLMAX .RESPTS),  Z(NLMAX  .RESPTS),  PI,  DELTHETA,  DEGTORAD 

REAL  THETA,C(5).R(5),H(5),DELZ(5),SHFT  .FACTOR, SQP.ZF 

REAL  ALPHA  1(5),  ALPHA(5),  B,  M,  RL, GAMMA 

COMPLEX  ER.UR 

CHARACTER*!  YN 

C  SET  CONSTANTS 
C 

PI=3. 1415927 
SQP=SQRT(PI) 

DEG  fORAD  =  PI/ 180. 

C 

C  DEFINE  START  AND  END  POINT  OF  EACH  RADOME  LAYER 
C 

Z(L,1)=0.0 

X(L,1)=R(L) 

Z(L,RHSPTS)=H(L) 

X(L,RESPTS)=0.0 

C . DEFINE  RADOME  LAYER  BY  VON  KARMAN  EQUATIONS  .... 

DELZ(L)=H(L)/RESOLV 

DO  33  I=2,RESPTS-1 

Z(L,0=Z(L,I- 1  )+DELZ(L) 

C  EQUATION  OF  RADOME  IS 
C  x=(  R/(SQRT(PI)) )  *  SQRT(  FZ  -  SIN(2*FZ)/2  ) 
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C  for  FZ  =  ACOS(  1  -  2  *  (H-z)/H  ) 

FZ  =  ACOS(  1 .0  -  2.0  *  ( H(L)-Z(L,I)  )  /  H(L) ) 
X(LJ)  =  R(L)/SQP  *  SQRT(  FZ  -  0.5*SIN(2*FZ) ) 

33  CONTINUE 


C . END  OF  RADOME  LAYER  .... 

RETURN 

END 

,****,*»**,**,,*.*****END  op  VQN  karman 
*********************************** 

C  POWER  SERIES  SHAPING  FUNCTION  SUBROUTINE 

SUBROUTINE  POWER(L,X,Z,H,R,F,RESOLV,RESPTS,NLMAX) 

INTEGER  IBIG,  I,  L,  SEPRATN,  NSHAPE.  NLMAX,  RESPTS,  RESOLV 
REAL  X (NLMAX .RESPTS),  Z(NLMAXJRESPTS),  PI.  DELTHETA,  DEGTORAD 
REAL  THETA,C(5)JR(5)JI(5)T)ELZ(5),SHFTEACTOR 
REAL  ALPHA1(5),  ALPHA(5),  B,  M,  RL, GAMMA,  F(NLMAX,7) 

COMPLEX  ER.UR 
CHARACTER*  1  YN 

C  SET  CONSTANTS 
C 

PI=3. 1415927 
DEGTORAD  =  PI/180. 

C 

C  DEFINE  START  AND  END  POINT  OF  EACH  RADOME  LAYER 
C 

Z(L,1)=0.0 

X(L,1)=R(L) 

Z(L,RESPTS)=H(L) 

X(L^ESPTS)=0.0 


C . DEFINE  RADOME  LAYER  BY  OGIVE  EQUATIONS  .... 
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DELZ(L)=H(L)/RES  OL  V 


C  Prompt  for  Shaping  factor  for  this  power  series  layer 
C 

WRITE(V) 

WRITE(*,43) 

43  FORMAT(8X,'  The  power  series  shaping  function  requires  a’, 
1/8X, ’shaping  factor,  F6,  to  be  used  in  the  exponential  in’, 
l/8X,'order  to  determine  the  radome  construction.', 

1//8X,'  When  F6  =  1.0  the  radome  shape  is  conical', 

1/8X,'  For  F6  =  0.75  the  radome  shape  approximates  the  ’, 
l/8X,7s!ewtonian  contour  (minimizes  drag).  Please  enter  F6:') 

READ(*,*)  F(L,6) 

DO  33  I=2,RESPTS-1 

Z(L,0=Z(L.I-1)+DELZ(L) 

C  EQUATION  OF  RADOME  IS  x=R  *  (  (H-z)/H  )**F6 
C 

X(L,I)  =  R(L)  *  ( (H(L)-Z(L,I))  /  H(L)  )**F(L,6) 

33  CONTINUE 


C END  OF  RADOME  LAYER  .... 

RETURN 

END 

**************** . END  OF  POWER  SERIES 

*•*••••••***•••••***•**•••••*••*••* 

C  PARABOLIC  SHAPING  FUNCTION  SUBROUTINE 

SUBROUTINE  PARACRV(L,X,Z,H,R,F,RESOLV .RESPTS  .NLMAX) 

INTEGER  IBIG,  I,  L,  SEPRATN,  NSHAPE,  NLMAX,  RESPTS,  RESOLV 
REAL  X (NLMAX .RESPTS),  Z(NLMAX,RESPTS),  PI,  DELTHETA,  DEGTORAD 
REAL  THETA,C(5)^(5),H(5),DELZ(5),SHFTf  ACTOR 
REAL  ALPHA  1(5),  ALPHA(5),  B,  M,  RL, GAMMA,  F(NLMAX,7) 

125 


COMPLEX  ER,UR 
CHARACTER*  1  YN 

C  SET  CONSTANTS 
C 

PI=3. 1415927 
DEGTORAD  =  PI/180. 

C 

C  DEFINE  START  AND  END  POINT  OF  EACH  RADOME  LAYER 
C 

Z(L,1)=0.0 

X(L,1)=R(L) 

Z(L,RESPTS)=H(L) 

X(L,RESPTS)=0.0 


C DEFINE  RADOME  LAYER  BY  PARABOLIC  EQUATION  .... 


C  Prompt  for  Shaping  factor  for  this  parabolic  layer 
C 

WRITE(*  ,*) 

WRITE/*  ,43) 

43  FORMAT(8X,'  The  parabolic  shaping  function  requires  a’, 
l/8X,'shaping  factor,  F7,  to  determine  the  curvature  of, 
l/8X,'the  parabolic  arc  for  the  radome  construction.' 

1//8X, 'Please  enter  F7:') 

READ(*,*)F(L,7) 

C . FOR  ROUNDED  PARABOLIC  RADOME  LAYER  .... 

C  ...  RADOME  IS  SCALED  TO  A  HEIGHT  =  SCALE  FUNCTION.C 

DELZ(L)=F(L,7)/RESOLV 


DO  33  1=2,RESPTS-1 


126 


Z(L,I)=Z(LJ-1)+DELZ(L) 


C  RECALL  EQUATION  OF  RADOME  IS  X*2=F7-Z,  X=SQRT(F7-Z) 
X(L,I>SQRT(  F(L.7)-Z(L4)  ) 

33  CONTINUE 

C 

C  SCALE  POINTS  TO  RADOME  DIMENSIONS 
C 

DO  44  I=2,RESPTS- 1 
X(L,I)=X(L,I)*R(L)/SQRT('  FC  ,,) ) 

Z(L,I)=Z(L,I)*H(L)/F<L,7' 

44  CONTINUE 

C . END  OF  ROUNDED  RADOME  LAYER  .... 


C . END  OF  RADOME  LAYER  .... 

RETURN 

END 

*************«**», ****END  0F  pARA  curve 
*********************************** 

C  POINTED  PARABOLIC  SHAPING  FUNCTION  SUBROUTINE 
SUBROUTINE  APEX(L,X,Z,H,R,F,RESOLV,RESPTS,NLMAX) 

INTEGER  IBIG,  I,  L,  SEPRATN,  NSHAPE,  NLMAX,  RESPTS,  RESOLV 
REAL  X(NLMAX,RESPTS),  Z(NLMAX^ESPTS),  PI,  DELTHETA,  DEGTORAD 
REAL  THETA, C(5)Jl(5),H(5),DELZ(5),SHFTJ;ACTOR 
REAL  ALPHA  1(5),  ALPHA(5),  B,  M,  RL, GAMMA,  F(NLMAX,7) 

COMPLEX  ER.UR 
CHARACTER*  1  YN 

C  SET  CONSTANTS 
C 
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PI=3. 1415927 
DEGTORAD  =  PI/180. 

C 

C  DEFINE  START  AND  END  POINT  OF  EACH  RADOME  LAYER 
C 

Z(L,1)=0.0 

X(L,1)=R(L) 

Z(L,RESPTS)=H(L) 

X(L,RESPTS)=0.0 

C . DEFINE  RADOME  LAYER  BY  POINTED  PARABOLIC  EQUATION  .... 


C  Prompt  for  Shaping  factor  for  this  parabolic  layer 
C 

WRITE(*  ,*) 

WRITE(*,43) 

43  FORMAT(8X,'  The  APEX  parabolic  shaping  function  requires  a’, 
l/8X,'shaping  factor,  F7,  to  determine  the  curvature  of, 
l/8X,'the  parabolic  arc  for  the  radome  construction  and  the', 
l/8X,’percent  of  the  radome  that  is  the  pointed  apex.', 
l//8X,'Please  enter  F7:‘) 

READ(*,*)F(L,7) 

C .  FOR  POINTED  RADOME  LAYER . 

C  ...  RADOME  IS  SCALED  TO  A  HEIGHT  =  MODIFIED  SCALE  FUNCTION E7+1 
DELZ(L)=(  F(L,7)+1.0  )/RESOLV 

D0  55I=2,RESPTS-1 

Z(L,I)=Z(L,I- 1  )+DELZ(L) 

C 

IF(  Z(L,I)  LE.(  F(L,7)  -1  ) )  THEN 
RECALL  EQUATION  OF  CONE  IS  X*2=F7-Z,  X=SQRT(F7-Z) 
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C 


X(L,I)=SQRT(  F(L,7)-Z(L,I) ) 

ELSE 

C  STRAIGHT  LINE  FORMULA  FOR  POINTED  END 
X(LJ)=(  ( F(L,7)  +1  )-Z(L,I)  )/2.0 
ENDIF 


55  CONTINUE 

C 

C  SCALE  POINTS  TO  RADOME  DIMENSIONS 
C 

DO  66  I=2,RESPTS-1 
X(L,I)=X(L,I)*R(L)/SQRT(  F(L,7) ) 
Z(L,I)=Z(L,I)* H(L)/(  F(L,7)+1.0 ) 

66  CONTINUE 

C . END  OF  POINTED  RADOME  LAYER  .... 


C . END  OF  RADOME  LAYER  .... 

RETURN 

END 

*****.*******, »»******END  QF  PARA  Ap£x 
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Appendix  C 


CORFU)  PROGRAM  CODE 


PROGRAM  CORFLD 

C  Program  Created  by  LT  RJ.  Vince  (Aug-Dec  89) 

C  Modifications  by  M.A.  Morgan  (Nov-Dec  89) 

C  Computing  core  E-field  at  each  point  along  a  50x50  radome  grid 
C  in  the  radome's  field  coordinate  system . 

C 

INTEGER  BELL,FF,NPTS,L.J.K,K1,NA,NINC 
INTEGER  MST ART .MSTOP.N  1  JVM  1  M.MFMD.M  1  M2, MM 
INTEGER  I  JIE4TM,KN,CA,CBNJN,NLP,SELECT,SELCT1 
REAL  ALPHA(  1 3),D  ALPHA(  1 3)  ANGINC(  1 3) 

REAL  THETA,PHipi,DTR,RMIN,RCORE 

REAL  XP,YP,ZP,R,COST,SINT,COSP,SINP,S,C,CMP,SMP,SA,SAMX 
REAL  R 1  ,R2,RSQ£M,TH,FI,RPTS  ,L  1 J 1  J>LOT(50,50),  A(3,3) 

REAL  P(40),DP(40),TP(6),PP(6) 

REAL  EOTM£OTE,E  I(8),EDTERR(50,50,2)1EDTER  1  ,EDTER2 
COMPLEX  EXP  1  £(8)  JEDOT (50.50^2)  ,EDOT  1  ,EDOT2 
COMPLEX  EP1,EP2,ET1£T2,ER1JER2JAY,RK,SRE,SR1J,KR 
COMPLEX  ERAD(50,50,2),ETH(50,50,2),EPH1(50,50,2) 

COMPLEX  ER(  1 6),UR(  16) ,COEF(70,80),CJ (40),DCJ (40) ,CPLOT 
CHARACTER*  1  YN 
CHARACTER*8  PRTDAT 
CHARACTER*  12  PRTDT 1  .PRTDT2  J>RTDT3 
CHARACTER*64  GRFLAB  ,HDR  1  ,HDR2,TITLE  1  ,TITLE2,DUM 
SLARGE 
C 

C  Setting  constants 
BELL=CHAR(7) 

FF  =  CHAR(12) 

PI  =  3.14159265359 
DTR=  PI/180. 

JAY=  (0.0, 1.0) 

C  Set  default  angle  of  rotation  of  planar  antenna  array  away  from  z  axis 

C  Set  default  aspect  ratio  for  plotting 
ASPECT  =  1.35 

C  Set  number  of  points  of  resolution  across  array 
NPTS=50 

RPTS=(  FLOAT(NPTS)  -  1.  )/2. 


C**lnitialize  arrays 
C  Initialize  the  E-field  matrices 
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on  on  no  on 


DO  1  L=1,NPTS 
DO  1  J-i  jNfPTS 
DC  1  K=l,2 
ERAD(J,L,K)=(0.,0.) 

ETH(JJLtK)=(0.,0.) 

EPHI(J,L,K)=(0.,0.) 

1  CONTINUE 

Prompt  for  input  data  filename  without  extension 

WRITE(*,14) 

READ(M01)  PRTDAT 
NCHAR=INDEX(PRTDAT,’ r)  -1 
PRTDT1  =  PRTDAT(l.NCHAR)  //  '.INT 
OPEN(  1 4  JFILE=PRTDT  1 ) 

TITLE  l=The  datafile  used  to  generate  fields  is  ’ 
TITLE  1  =TITLE  1  ( 1 :40)  //  PRTDT 1 

PRTDT2=  PRTDAT(1:NCHAR)  //  '.DAT' 
OPEN(  1 5  JTLE=PRTDT2) 


Write  header  from  structure  data  file 

READ(14,101)  HDR1 
READ(14,101)  HDR2 
READ(14,101)  GRFLAB 
WRITE(15t101)  HDR1 
WRITE(15,101)  HDR2 
WRITE(15, 101)  GRFLAB 
WRITER,  101  )HDR1 
WRITER,  101 )  HDR2 
WRITE(*, 101)  GRFLAB 


Read  interior  core  Er.Ur 

READ(14,507)  ER(1),DUM 
READ(  14,507)  UR(1),DUM 
SRE=CSQRT(ER(1)) 

SRU=CSQRT(UR(1)) 

KR=SRE*SRU 
WRITER,*) 

WRITE(*,*)'SRE=  ',SRE 
WRITE(*,*)’SRU=  ',SRU 
WRITE(*,*)’KR  =  ',KR 

WRITE(15  *)'  SRE=',SRE,’;  SRU=',SRU,’;  KR  =  ,KR 

C  Read  minimum  diameter  of  missile  radome  core  expansion 
C  Assuming  in  Ko*r  Wavenumber  Units 
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non  on  o  non  on  on  on 


c 

RE AD(  14,503)  RMIN.DUM 
RCORE=RMIN 


Read  number  of  incident  field  angles  (.LE.  5) 

READ(  14,500)  NA,DUM 
NDMC=2*NA 


Read  incident  angles 

IPZ=0 

IMZ=0 

DO  11 1=1  ^NA 

READ(14,5 11)  DANGLE 

Checking  for  +/-  Z-Axis  Incident  Fields 

BF(DANGLE.EQ.O.O)  IPZ=I 
1F(DANGLE.EQ.  180.0)  IMZ=I 

DALPHA(I)=  DANGLE 

EMRAD  uses  Poynting  vector  angle  with  Z-axis  reference 
ANGle  of  incidence  =  180  -  D ALPHA  from  EMRAD 

ANGINC(I)=  180-  DANGLE 


WRITE(*,*) '  Incoming  Inc.  Angle  #',I 

WRITE(*,*) '  EMRAD  Poynting  angle  of  incidence  =  ,DALPHA(I) 

WRITE(*,*) '  Radome  incident  angle  =  ,ANGINC(1) 

WRITE(*,*) 

WRITE(15,*)  ’  Incoming  Inc.  Angle  =  ’.DALPHA(l) 

Set-up  array  of  angles  in  radian  measure 

ALPHA(I)=DTR*DALPHA(I) 

1 1  CONTINUE 
NAI=1 

IF  (NA.EQ.1)  GO  TO  13 
Determine  angle  of  interest 

12  CONTINUE 

WRITE(*,*) '  Select  Angle  of  Interest  By  Number  ’ 

WRITE(*,*) 

READ(*,*)  NAI 

IF  ((NAI.GE.l).OR.(NAI.LE.NA))  GO  TO  13 
WRITE(*,*) '  That  Number  is  Not  Available  ...  Try  Again’ 
WRITER,*) 
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no  no  no 


GOTO  12 
13  CONTINUE 


Checking  if  Selected  Field  is  +/-  Z-Axis  Incident 
IM1=0 

IF((NAI.EQ.IPZ)OR.(NAI£Q.IMZ))  IM1=1 
WRITER, *) 

WRITE(*,*) '  The  angle  of  consideration  is  angle  #’,NAI 
1  ,’  =  ',ANGINC(NAI),'  degrees' 

WRITE(V) 

WRITE(15,*) '  The  angle  of  consideration  is  angle  #’,NAI 
1  ,'  = ',ANGINC(NAI),' degrees' 

WRITE(15,*) 

Transformation  of  coordinates  to  gimble  mounted  planar  array 

WRITE(*,*)'  The  planar  array  is  set  by  default  to  lie  in  the’ 
WRITE(*,*)'  x-y  plane,  with  boresight  in  the  z-direction.' 
WRITER,*) 

WRITE(*,*)'  Do  you  want  to  change  the  array  boresight  ?’ 
WRITER,*)'  ( "Y"  or  "N" ) ' 

IFLG=0 
DTH=0.0 
READ(M01)  YN 

IF(YN.EQ.T  .OR.YN.EQ.'y')  THEN 
WRITE(V) 

WRITE(*,*)  "Enter  Boresight  Theta  Angle  (in  degrees)’ 
READ(V)  DTH 
TH  =  DTH  *  DTR 
WRITE(V) 

WRITE(*,*)  'Enter  Boresight  Phi  Angle  (in  degrees)’ 
READ(*,*)  DR 
FI  =  DFI  *  DTR 

IFLG=1 

Loading  Transformation  Matrix  from  Array  to  Core  Coordinates 

CT=COS(TH) 

ST=SIN(TH) 

CP=COS(FI) 

sp=siN(n) 

A(1,1)=CT*CP**2.0+SP**2.0 
A(  1 ,2)=CP*  SP*(CT  - 1 .0) 

A(1,3)=ST*CP 

A(2, 1  )=CP*  S  P*  (CT  -1.0) 

A(2,2)=CT‘SP**2.(H-CP**2.0 

A(2,3)=ST*SP 

A(3,1)=-ST*CP 
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on  on 


A(3,2)=-ST*SP 

A(3,3)=CT 

END  IF 


Get  modal  information 

RE AD(  14,505)  MST ART .MSTOPDUM 

WRITE(* ,*)’  MST ART=',MSTART,'  MSTOP=',MSTOP 

WRITE(15,*) '  MSTART=',MSTART,'  MSTOP=',MSTOP 

Loop  modals  to  read  modal  coefficients 

WRITE!*,*)  ’Reading  Coefficients  and  Constructing' 

WRITE!*,*)  "Fields  for  Each  Azimuthal  "m"  Mode' 

WRITE!*,*) 

IF((MSTART.EQ.0).AND.(IM1.EQ.1))  THEN 

WRITE(*,*)  Terminating  M-Loop  at  m=l  for  +/-  Z-Axis  Incidence' 

WRITE!*,*) 

MSTOP=l 

ENDIF 

DO  777  MD=MSTART,MSTOP 

READO4.500)  M.DUM 
WRITE!  15,*)'  M=  ',M 
DF(M£Q.MD)  GO  TO  220 

WRITE!*,*)  'Error  in  Reading  "m"  -  Check  Data  File' 

STOP 

220  CONTINUE 
WRITE(V) 

WRTTE(*,*)  ’Azimuthal  Mode  "m"  =  ',M 
WRITE(**) 

BF(M.EQ.0)THEN 
EM  =  1.0 
MF  =  1 
ELSE 
EM  =  2.0 
MF  =  0 
ENDIF 

WRITE(*,*)  ’Reading  Coefficients' 

READ!  14,500)  NMl.DUM 
DO  665  ID=I,NA 
READ!  14,500)  I DUM 
IF(ID.EQ.I)  GO  TO  221 

WRITE!*,*)  'Inc  Angle  Index  Read  Error  -  Check  Data  File’ 

STOP 

221  CONTINUE 
WRITE(15*)'  1=  ’,1 
ITM=I 


134 


noon  n  n 


DO  222  KN=1,NM1 
CA=KN 
CB=CA+NM1 

READ(14,555)  COEF(CA,ITM),COEF(CB,ITM),DUM 
WRTTE(15,*)  KN,COEF(CAJTM),COEF(CBJTM)X)UM 
222  CONTINUE 
ITE=I+NA 
DO  333  KN=lJsTMl 
CA=KN 
CB=CA+NM1 

READ(14,555)  COEF(CA,ITE),COEF(CBJTE)>DUM 
WRITE(15,*)  KN,COEF(C A JTE) ,COEF(CB  ,ITE) ,DUM 
333  CONTINUE 
665  CONTINUE 

C  Skipping  Field  Generation  for  m=0  if  +/-  Z- Axis  Incidence 
IF((M.EQ.0).AND.(IM1.EQ.1))  GO  TO  777 

WRITE(*,*)  'Generate  Fields  Across  Array' 

DX=RCORE/RPTS 

DY=DX 


DO  55  J=1NPTS 
DO  33  L=1NPTS 


Determine  local  planar  array  coordinates 

Jl=FLOAT(J-l)-RPTS 

LI=FLOAT(L-l)-RPTS 

XP=J1*DX 

YP=L1*DY 

ZP=0.0 

If  Array  Tilted  Then  Transform  Coordinates 
From  local  Planar  array  coordinates 

to  global  Core  coordinates.... 

IF(IFLG.EQ.l)  THEN 
XC=A(  1 ,1  )*XP+A(  1 ,2)*YP+A(  1 ,3)‘ZP 
YC=A(2,1)*XP+A(2,2)*YP+A(2,3)*ZP 
ZC=A(3,l)*XP+A(3r2)*YP+A(3,3)*ZP 
ELSE 
XC=XP 
YC=YP 
ZC=ZP 
END  IF 

Determine  if  R  <=  RCORE,  Else  Skip  Field  Computation 

RSQ=XC*XC+YC* YC+ZC*ZC 
R=SQRT(RSQ) 
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no  oo  on  no  no  on 


IF  (R.GT.RCORE)  GO  TO  32 

Determine  spherical  Bessel  function  scaled  radius  argument 
RK=KR*R 

Determine  Complex  Spherical  Bessel  function  and  Derivative 

NJN=M+NM  1 +MF 

CALL  CSBSL(RK,N  JN,CJ  JXU) 

Compute  Core  Coordinates  and  Functions 

THETA=ACOS(ZC/R) 

PHI=ARCTAN2(YC,XQ 
C=COS(THET  A) 

S=SEN(THETA) 

CMP=EM*COS(M*PHI) 

SMP=EM*SIN(M*PHI) 

Determine  Legendre  Polynomial  and  Derivative 

NLP=NM1+MF 

CALL  LPAD(M,NLP,C,P,DP) 

Consider  only  the  incident  angle  of  interest 


I=NAI 

ITM=I 

ITE=ITM+NA 

Initialize  additive  KN-sum  fields  to  zero 


EP1=(0.,0.) 

EP2=(0.,0.) 

ET1=(0.,0.) 

ET2=(0.,0.) 

ER1=(0.,0.) 

ER2=(0.,0.) 


DO  22  KN=1,NM1 
N=KN+M+MF-1 
NJ=N+1 
NL=KN+MF 
CA=KN 
CB=CA+NM1 

ET1=ET1+C0EF(CA,ITM)*M*CJ(NJ)*P(NL) 
1  +COEF(CB,ITM)*DCJ(NJ)*DP(NL) 

EP 1  =EP  1  +COEF(CA,ITM)*CJ(NJ)*  DP(NL) 


1 


+COEF(CBJTM)*M*DCJ(NJ)*P(NL) 

ER 1  =ER  1  +COEF(CB  JTM)*N*(N+1)*CJ(NJ)*P(NL) 


ET2=ET2+COEF(CA,ITE)*M*CJ(NJ)*P(NL) 

1  +COEF(CBJTE)*DCJ(NJ)*DP(NL) 

EP2=EP2+COEF(CAJTE)*CJ(NJ)*DP(NL) 

1  +COEF(CB  JTE)*  M*  DCJ  (NJ)*  P(NL) 

ER2=ER2+COEF(CB4TE)*N*(N+l)*a(NJ)*P(NL) 

22  CONTINUE 


ETH(J  ,L,1  )=ETH(J,L,  1  )-J  AY*SRU*CMP*ET1/R 
EPHI(J4-,1)=EPH1(JJL,1>+JAY*SRU*SMP*EP1/R 
ERAEKJX,  1  )=HR  AEKJJ-.  1  )-J  A  Y  •  S  *CMP*ER  1 /(SRE*RSQ) 

ETH(  J  JL^)=ETH(I  ,L,2>+  SRU*  S  MP*  ET2/R 

EPHI(JU,2)=EPHI(JJL^)+SRU*CMP*EP2/R 

ER  AD(J  ,L,2)=ER  AD(  J  J-,2)+  S  *  S  MP*  ER2/(S  RE  *  R  S  Q) 

32  CONTINUE 

33  CONTINUE 

* 

55  CONTINUE 
C 

C—  End  of  loops  across  antenna  plane  -  ■= 

111  CONTINUE 
C 

C**  End  of  modal  loop  for  field  analysis 


Q+ ************** 

C  Determine  magnitude  and  phase  of  Computed  field  in  the  direction 
C  of  the  theoretical  field  vector,  and  the  magnitude  of  the  error. 

C 


SINA=SIN(ALPHA(NAI)) 
COS  A=COS(ALPHA(NAI)) 


C  Loop  through  array  to  determine  field  components 
C 

DO  634  J=1  J'TPTS 
DO  633  L=1  JnTTS 

C  Determine  local  planar  array  coordinates 
C 

Jl=FLOAT(J-l)-RPTS 

Ll=FLOAT(L-l)-RPTS 

XP=J1*DX 

YP=L1*DY 

ZP=0.0 
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C  If  Array  Tilled  Then  Transform  Coordinates 
C  From  Array  to  Core  .... 

C 

IF{IFLG£Q.l)  THEN 
XC= A(  1,1)*  XP+A(  1 ,2)*YP+A(  1 3)*ZP 
YC=A(2,1)*XP+A(2,2)*YP+A(23)*ZP 
ZC=A(3,1)*XP+A(33)*YP+A(33)*Z 
ELSE 
XC=XP 
YC=YP 
ZC=ZP 
END  IF 

RSQ=XC*  XC+ Y  C*  YC+ZC*  ZC 
R=SQRT(RSQ) 

C  Determine  if  R  <=  RCORE,  Else  Zero  Fields 
C 

IF  (R.LE.RCORE)  THEN 

C  Compute  Core  Coordinates  and  Functions 
C 

THETA=ACOS(ZC/R) 

PH1=ARCTAN2(YC,XC) 

CT=COS(THETA) 

ST=SIN(THETA) 

CP=COS(PHI) 

SP=SIN(PHI) 

C  NOTE:  R  IS  GIVEN  IN  TERMS  OF  KO*R 
C  =  KO  IS  INCLUDED  IN  R,  DEFINE  KO=  1 .0 
C 

KO=1.0 

C  TM  and  TE  Plane  waves  are  of  unit  magnitude 
C 

EOTE=  1 .0 
EOTM=  1 .0 


C  Complex  incident  field  is  exponential(-j*ko*p.r) 

C  use  field  direction  for  incident  angle  1 80  to  90  degrees 
C 

EXP1  =  CEXP(  -JAY*KO*(  XC*SINA+ZC*COSA  ) ) 

C  Compute  plane  wave  field  unit  vectors  and  components 
C 

C  E  TM  =  cos(alpha)Ux  -  sin(alpha)Uz 
C 

C  E  TM  T 
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o  o 


EI(1)  =  COSA*CT*CP  -  SINA*(-ST) 

E(1)  =  E0TM*EXP1*EI(1) 

C  ETMP 

EI(2)  =  COSA*(-SP) 

E(2)  =  E0TM*EXP1  *EI(2) 

C  ETMR 

EI(3)  =  COSA*ST*CP  -  SENA*(CT) 

E(3)  =  E0TM*EXP1*EI(3) 

C  E  TM  TOT 

E(4)=SQRT(CABS(E(1))**2.0+CABS(E(2))**2.0+CABS(E(3))**2.0) 
E  TE  =  Uy 


C  ETET 

EI(5)  =  CT*SP 

E(5)  =  EOTE*EXPl*EI(5) 

C  ETEP 

EI(6)  =  CP 

E(6)  =  EOTE*EXPl*EI(6) 

C  ETER 

EI(7)  =  ST*SP 

E(7)  =  EOTE*EXPl*EI(7) 

C  E  TE  TOT 

E(8)=SQRT(CABS(E(5))**2.0+CABS(E(6))**2.0+CABS(E(7))**2.0) 

*»***»»**» 


C  Find  magnitude  and  phase  of  Computed  field  in  the  direction 
C  of  the  theoretical  field  vector 

C  E  .  ei  /  Ei  =  (  eth*ETH+ephi*EPHl+erad*ERAD  )/EO*EXPl 
C 

EDOTl=EI(l)*ETH(JJL,l)+EI(2)*EPHI(J^,l)+EI(3)*ERAD(J(L>l) 
EDOT(J  JL,  1  )=EDOTl/(EOTM*EXP  1 ) 

EDOT2=EI(5)*ETH(JJLv2)+EI(6)*EPHI(JT.,2)+EI(7)*ERAD(J^,2) 
EDOT(J  ,L,2)=EDOT2/(EOTE*EXP  1 ) 

C  Find  rel.  magnitude  of  Computed  field  perpendicular  to  the  direction 
C  of  the  theoretical  field  vector,  the  error  component  in  the 
C  Computed  field. 

C 

C  Rel  error  =  I E  -  (E.ei)  ei  I  /  EO 
C 

C  =  sqrt{  Eth  -  (E.ei)*ethl**2  +  lErad  -  (E.ei)*eradl**2 
C  +  Ephi  -  (E.ei)*ephil**2  }/  EO 
C 

C  Field  error  magnitude  is  computed 
C 


EDTER1  =  CABS(ETH(J,L,1)-EDOT1*EI(1))**2.0 
1  +  CABS(EPHI(JX,1)-EDOT1*EI(2))**2.0 

1  +  CABS(ERAD(J,L,1)-EDOT1*EI(3))**2.0 
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EDTERR(JJ-,  1  )=SQRT(EDTER  1  )/(EOTM) 

EDTER2=  C ABS(ETH(J  J-,2)  -EDOT2*EI(5))**2.0 
1  +  CABS(EPHI(J,L,2)-EDOT2*EI(6))**2.0 

1  +  CABS(ERAD(J,L,2)-EDOT2*EI(7))**2.0 

EDTERR(J  ,L,2)=SQRT  (EDTER2)/(E0TE) 


ELSE 


C  R  >  RCORE  ,  Set  field  components  to  zero 
C 

EDOT(J,L,1)=0.0 
EDOT(J  ,L,2)=0.0 
EDTERR(JJL,1)=0.0 
EDTERR(J  J-,2)=0 .0 
ENDIF 

633  CONTINUE 

634  CONTINUE 

£********** 


WRITE(*  *)  'SPHERICAL  FIELDS  COMPUTED' 

WRITE(V) 

C  Optional  Listing  of  Coordinates  and  E-Fields 
C 

WRITEf*,*)  'Check  Values  of  Coordinates  and  E-Field  ?  (Y/N):' 
READ(M01)  YN 

IF((YN.EQ.T')OR.(YN.EQ.'y'))  THEN 
757  CONTINUE 

WRITE(*,*)  'Enter  Array  Triple  Index  to  Check  (J,L,I)' 
WRITE(* ,*)  ’  J=  x-point  1,2,..,\NPTS 

WRITE(**)'  L=  y-point  1,2,..,’,NPTS 

WRITE(*  *) '  1=  1  for  TM,  2  for  TE' 

WRITE(*  ,*) 

READ<*  ,*)JL.I 
EPM=CABS(EPHI(JL,D) 

ETM=CABS(ETH(J  ,L,I)) 

ERM=CABS(ERAD(J  ,L  J)) 

ETOT»SQRT(EPM*EPM+ETM*ETM+ERM*ERM) 


C  Determine  local  planar  array  coordinates 
C 

Jl=FLOAT(J-l)-RPTS 
Ll=FLOAT(L-l)-RPTS 
XP=J1*DX 
YP=L1*DY 
ZP=0.0 

1 

C  If  Array  Tilted  Then  Transform  Coordinates 
C  From  Array  to  Core  .... 
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c 

IF(IFLG£Q.l)  THEN 
XC= A(  1 , 1  )*XP+ A(  U)*  YP+ A(  1 ,3)*ZP 
YC=A(2,1)*XP+A(2,2)*YP+A(2,3)*ZP 
ZC=A(3,1)*XP+A(3,2)*YP+A(3,3)*ZP 
ELSE 
XC=XP 
YC=YP 
ZC=ZP 
ENDIF 

WRITE(*,*)  'XP,YP.ZP:',XP,YP,ZP 
WRITER,*)  'XC,YC,ZC:'.XC,YC,ZC 
WRITE(**) 

WRITER, *)  ’IE-RADI  ='3RM 
WRITE(V) 'E-THI  ='ETM 
WRITE(*,*)  'E-PHII  ='£PM 
WRITE(*  ,*)  'IE-TOTI  ='£TOT 
WRITE(*,*) 

WRITE(*,*)  ’Look  at  Another  Point  ?  (Y/N):’ 

REACH*, 101)  YN 

IF((YN.EQ.'Y’).OR.(YN.EQ.,y'))  go  to  757 
ENDIF 

£  ********************************************************************** 
C  Field  representation  selection  menu 
^  *********************************** 

c 

700  WRITE(*,701) 

Q  ****************************************************************** 

701  f orm at(////////////5 X ,'COR FL D  OUTPUT  GRAPHICS  PRESENTATION V8X, 
l'CORE  FIELD  REPRESENTATION  SELECTION  MENUV/5X, 

2Please  select  the  type  of  field  representation  of  interest 
375X, 

4'CORE  will  display  the  following  field  representations:7/8X, 

5'1.  TM  INCIDENT  FIELD,  E-THETA78X, 

6'2.  TM  INCIDENT  FIELD,  E-PHI  78X, 

7'3.  TM  INCIDENT  FIELD,  E-RADIAL  78X, 

84.  TM  INCIDENT  FIELD,  E-TOTAL  FIELD  MAGNITUDE  7/8X, 

9'5.  TE  INCIDENT  FIELD,  E-THETA',/8X, 

1  '6.  TE  INCIDENT  FIELD,  E-PHI  78X, 

27.  TE  INCIDENT  FIELD,  E-RADIAL  '^X, 

3'8.  TE  INCIDENT  FIELD,  E-TOTAL  FIELD  MAGNITUDE  7/8X, 

4-9.  COMPARE  COMPUTED  FIELD  TO  INCIDENT  FIELD7/8X, 

510.  CHANGE  ASPECT  RATI07/8X, 

6  0.  FINISHED  WITH  THIS  ANGLE  OF  INCIDENCE7/5X, 

6'Indicate  your  selection  entering  ”l","2",...,"10",or  ”0" 

7',5X) 

p  ********************************** 


READ(V)  SELECT 
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C  Check  response  to  menu 


IF  (SELECTEQ.l)  THEN 

Q  ****************************************************************** 

C  DISPLAY  TM  INCIDENT  FIELD,  E-THETA 

C  CHECK  FOR  DESIRED  REPRESENTATION  OF  FIELD,  MAGNITUDE/PHASE 
WRITE(*,715) 

Q  ******** ************ *•••••*•*******•»***••******** ****** ********** 

715  format(////////////5X,'CORE  FIELD  PRESENT ATIONS75X, 
l'You  have  selected  the  following  field  elemem7/5X, 

1 1 .  TM  INCIDENT  FIELD,  E-THETA7/8X, 

3This  field  can  be  plotted  in  PHASE  or  MAGNITUDE,’ ,//5X, 

4Please  select  the  type  of  field  representation  of  interest 
578X, 

6'1.  PHASE78X, 

7'2.  MAGNITUDE'ySX, 

8'Indicate  your  selection  entering  "1"  or  "2" 

97/5X) 

p  ********************************** 

READ(*,*)  SELCT1 

C  ASSIGN  AS  INDICATED  BY  RESPONSE  TO  MENU 

DO  712  J=1,NPTS 
DO  71 1  L=1 2MPTS 
CPLOT=ETH(JL,l) 

IF  (SELCTLEQ.l)  THEN 

PLOT(J^)=(ARCTAN2(IMAG(CPLOT),REAL(CPLOT)))/DTR 

ELSE 

PLOT(J  JL)=C  ABS(CPLOT) 

END  IF 

711  CONTINUE 

712  CONTINUE 
C 

IF  (SELCTLEQ.l)  THEN 

TITLE2  =  ' TM  INCIDENT  FIELD,  E-THETA  PHASE’ 

ELSE 

T1TLE2  =  ’  TM  INCIDENT  FIELD,  E-THETA  MAGNITUDE’ 

ENDIF 

C 

p  ****************************************************************** 


ELSEIF  (SELECT.EQ.2)  THEN 

Q  ****************************************************************** 

C  DISPLAY  TM  INCIDENT  FIELD,  E-PHI 

C  CHECK  FOR  DESIRED  REPRESENTATION  OF  FIELD,  MAGNITUDE/PHASE 
WRITER, 725) 

p  ****************************************************************** 
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725  format(////////////5X,’CORE  FIELD  PRESENTATIONS',^, 
lTou  have  selected  the  following  field  element' .//5X, 

2'2.  TM  INCIDENT  FIELD,  E-PM7/8X, 

3 This  field  can  be  plotted  in  PHASE  or  MAGN1TUDE,7/5X, 

4 Please  select  the  type  of  field  representation  of  interest 
5'JSX, 

61.  PHASE78X, 

7'2.  MAGNITUDE' y8X , 

8'Indicate  your  selection  entering  Tor  ”2" 

97/5X) 

Q  ********************************** 

RE  ADC*,*)  SELCT1 

C  ASSIGN  AS  INDICATED  BY  RESPONSE  TO  MENU 

DO  722  J=1,NPTS 
DO  721  L=1NPTS 
CPLOT=EPHI(J,L,l) 

IF  (SELCT1.EQ.1)  THEN 

PLOT(JJL)=(ARCTAN2(IMAG(CPLOT),REAL(CPLOT)))/DTR 

ELSE 

PLOT(J,L)=CABS(CPLOT) 

ENDIF 

721  CONTINUE 

722  CONTINUE 
C 

IF  (SELCT1.EQ.1)  THEN 

TTTLE2  =  ’  TM  INCIDENT  FIELD,  E-PHI  PHASE  ’ 
ELSE 

TTTLE2  =  '  TM  INCIDENT  FIELD,  E-PHI  MAGNITUDE' 
ENDIF 


Q  ****************************************************************** 

ELSEIF  (SELECT.EQ.3)  THEN 

Q  ****************************************************************** 

C  DISPLAY  TM  INCIDENT  FIELD,  E-RADI AL 

C  CHECK  FOR  DESIRED  REPRESENTATION  OF  FIELD,  MAGNITUDE/PHASE 
WRITE(*,735) 

Q  ****************************************************************** 

735  form at(////////////5 X ,’CORE  FIELD  PRESENT ATIONS75X, 
l’You  have  selected  the  following  field  eiement7/5X, 

2’3.  TM  INCIDENT  FIELD,  E-RADI  AL7/8X, 

3This  field  can  be  plotted  in  PHASE  or  MAGNTTUDE,7/5X, 

4Please  select  the  type  of  field  representation  of  interest 
578X, 

61.  PHASE78X, 

72.  MAGNITUDE78X , 

8  Indicate  your  selection  entering  "1"  or  "2" 


c 


97/5X) 

********************************** 

RE  ADC,*)  SELCT1 

C  ASSIGN  AS  INDICATED  BY  RESPONSE  TO  MENU 

DO  732  J=1,NPTS 
DO  731  L=1JVPTS 
CPLOT=ERAD(J,L.l) 

IF  (SELCTl£Q.l)  THEN 

PLOT(JEMARCTAN2(IMAG(CPLOT),REAL(CPLOT)))/DTR 

ELSE 

PLOT ( J  JL)=C  AB  S  (CPLOT) 

ENDIF 

731  CONTINUE 

732  CONTINUE 
C 

IF  (SELCT1.EQ.1)  THEN 

TTTLE2  =  ’  TM  INCIDENT  FIELD,  E-RADIAL  PHASE  ' 

ELSE 

TITLE2  =  •  TM  INCIDENT  FIELD,  E-RADIAL  MAGNITUDE1 
ENDIF 

C 

Q  ****************************************************************** 

ELSEIF  (SELECT.EQ.4)  THEN 

Q  •**••**•*•*•*•*•*•••**••**•*•*******»**•*•***••****••*•*******•*** 

C  DISPLAY  TM  INCIDENT  FIELD,  E-TOTAL  FIELD  MAGNITUDE 

C  ASSIGN  AS  INDICATED  BY  RESPONSE  TO  MENU 

DO  742  J=1,NPTS 
DO  74 1  L=1,NPTS 
EPM=CABS(EPHI(J  JL,  1 )) 

ETM=C  ABS  (ETH(J  E,  1 )) 

ERM=CABS(ERAD<J,L,1)) 

PLOT(JE)=SQRT(EPM*EPM+ETM*ETM+ERM*ERM) 

741  CONTINUE 

742  CONTINUE 

TITLE2  =  1  TM  INCIDENT  FIELD,  E-TOTAL  FIELD  MAGNITUDE  ’ 

Q  ****************************************************************** 

ELSEIF  (SELECT.EQ.5)  THEN 

Q  ****************************************************************** 

C  DISPLAY  TE  INCIDENT  FIELD,  E-THETA 

C  CHECK  FOR  DESIRED  REPRESENTATION  OF  FIELD,  MAGNITUDE/PHASE 
WRITE(*,755) 

( '  ****************************************************************** 
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755  form ax(////////////5 X ,'CORE  FIELD  PRESENTATIONS75X, 
l'You  have  selected  the  following  field  element7/5X, 

2’5.  TE  INCIDENT  FIELD,  E-THETA7/8X, 

3This  field  can  be  plotted  in  PHASE  or  MAGNITUDE, 7/5X, 
4 Please  select  the  type  of  field  representation  of  interest 
578X, 

61.  PHASE78X, 

77.  MAGNITUDE’y8X , 

8'Indicate  your  selection  entering  "1"  or  "2" 

97/5X) 

"»  ********************************** 


RE  ADC,*)  SELCT1 

C  ASSIGN  AS  INDICATED  BY  RESPONSE  TO  MENU 

DO  752  J=1,NPTS 
DO  751  L=1,NPTS 
CPLOT=ETH(JL,2) 

IF  (SELCT1.EQ.1)  THEN 

PLOT(JJL)=(ARCTAN2(IMAG(CPLOT),REAL(CPLOT)))/DTR 

ELSE 

PLOT(J  ,L)=C  AB  S(CPLOT) 

ENDIF 

751  CONTINUE 

752  CONTINUE 
C 

IF  (SELCT1.EQ.1)  THEN 

TITLE2  =  ' TE  INCIDENT  FIELD,  E-THETA  PHASE’ 

ELSE 

TITLE2  =  '  TE  INCIDENT  HELD,  E-THETA  MAGNITUDE' 
ENDIF 


Q  ****************************************************************** 

ELSEIF  (SELECT.EQ.6)  THEN 

£  ****************************************************************** 

C  DISPLAY  TE  INCIDENT  FIELD,  E-PHI 

C  CHECK  FOR  DESIRED  REPRESENTATION  OF  FIELD,  MAGNITUDE/PHASE 
WRITE(*,765) 

Q  ****************************************************************** 

765  f orm a t (//////////// 5 X ,’COR E  FIELD  PRESENTATIONS'.^, 
l’You  have  selected  the  following  field  element7/5X, 

2'6.  TE  INCIDENT  FIELD,  E-PHI7/8X, 

3This  field  can  be  plotted  in  PHASE  or  MAGNITUDE,7/5X, 

4Please  select  the  type  of  field  representation  of  interest 
578X, 

61.  PHASE78X, 

77.  MAGNITUDE’ ,/8X, 

8'Indicate  your  selection  entering  "1"  or  "2" 
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97/5X) 

Q  ********************************** 

READ{*  *)  SELCT1 

C  ASSIGN  AS  INDICATED  BY  RESPONSE  TO  MENU 

DO  762  J=1  NPTS 
DO  761  L=1,NPTS 
CPLOT=EPHI(JU,2) 

IF  (SELCTl£Q.l)  THEN 

PLOT(J,L)=(ARCTAN2(IMAG(CPLOT),REAL(CPLOT)))/DTR 

ELSE 

PLOT ( J  JL)=C  AB  S(CPLOT) 

ENDIF 

761  CONTINUE 

762  CONTINUE 
C 

IF  (SELCTLEQ.l)  THEN 

TTTLE2  =  '  TE  INCIDENT  FIELD,  E-PHI  PHASE  ' 

ELSE 

TTTLE2  =  '  TE  INCIDENT  FIELD,  E-PHI  MAGNITUDE  ' 

ENDIF 

C 

Q  •*•••*••**•***•••**•***•••*•*•••***••••*•••***»*•••**•*•*••••***•* 

ELSEIF  (SELECTEQ.7)  THEN 

q  ****************************************************************** 

C  DISPLAY  TE  INCIDENT  FIELD,  E-RADIAL 

C  CHECK  FOR  DESIRED  REPRESENTATION  OF  FIELD,  MAGNITUDE/PHASE 
WRITE(*,775) 

Q  ****************************************************************** 

775  format(////////////5X;CORE  FIELD  PRESENT AT10NS75X, 
l'You  have  selected  the  following  field  element7/5X, 

27.  TE  INCIDENT  FIELD,  E-RADIAL7/8X, 

3 This  field  can  be  plotted  in  PHASE  or  MAGNITUDE.7/5X, 

4 Please  select  the  type  of  field  representation  of  interest 
578X, 

67.  PHASE'^X, 

77.  MAGNITUDE78X, 

8'Indicate  your  selection  entering  "1"  or  "2" 

97/5X) 

Q  ********************************** 

READ(*  ,*)  SELCT1 

C  ASSIGN  AS  INDICATED  BY  RESPONSE  TO  MENU 

DO  772  J=1,NPTS 
DO  771  L=1NPTS 
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CPL0T=ERAD(J,L,2) 

IF  (SELCT1.EQ.1)  THEN 

PL0T(JJL)=(ARCTAN2(IMAG(CPL0T),REAL(CPL0T)))/DTR 

ELSE 

PLOT(J,L)=CABS(CPLOT) 

END  IF 

771  CONTINUE 

772  CONTINUE 
C 

IF  (SELCTl.EQ.l)THEN 

TITLE2  =  ' TE  INCIDENT  FIELD,  E-RADIAL  PHASE’ 

ELSE 

TTFLE2  =  ’  TE  INCIDENT  FIELD,  E-RADIAL  MAGNITUDE’ 

ENDIF 

C 

p  ****************************************************************** 


ELSEIF  (SELECT.EQ.8)  THEN 

£  ****************************************************************** 

C  DISPLAY  TE  INCIDENT  FIELD,  E-TOTAL  FIELD  MAGNITUDE 

DO  782  J=1,NPTS 
DO  781  L=1,NPTS 
EPM=CABS(EPHI(J,L,2)) 

ETM=CABS(ETH(JU,2)) 

ERM=CABS(ERAD(J,L,2)) 

PLOT(J,L)=SQRT(EPM*EPM+ETM*ETM+ERM*ERM) 

781  CONTINUE 

782  CONTINUE 

TITLE2  =  ’ TE  INCIDENT  FIELD,  E-TOTAL  FIELD  MAGNITUDE  ' 

****************************************************************** 

ELSEIF  (SELECT.EQ.9)  THEN 

a***************************** *»**••**•* *•••***•**•*•***•***•••**••*•* 


Incident  Field  comparison  representation  selection  sub  menu 

****************************** ***************** ************* 


800  WRITE(*,801) 

-  ****************************************************************** 

801  f orm at (//////////// 5 X ,’CORE  COMPARISON  TO  IDEAL  FIELD7/5X, 

1  'CORE  FIELD  REPRESENTATION  SELECTION  MENU  27/5X, 

2Please  select  the  type  of  field  representation  of  interest 

375X, 

4'CORE  will  display  the  following  field  representations:7/8X, 

5T .  TM  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field78X, 

6  2.  TM  INCIDENT  FIELD,  Scaled  Error  component  7//8X, 

73.  TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  FicldTSX, 

8'4.  TE  INCIDENT  FIELD,  Scaled  Error  component  7//8X, 

9'Indicate  your  selection  entering  "  1  ",”2”,"3",  or  ’’4" 
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17/5X) 

Q  ********************************** 

READ(*  ,*)  SELCT2 

C  Check  response  to  menu 

****** 


IF  (SELCT2£Q.l)  THEN 

Q  ****************************************************************** 

C  DISPLAY  TM  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field 
C  CHECK  FOR  DESIRED  REPRESENTATION  OF  FIELD,  MAGNITUDE/PHASE 
WRITE(*,815) 

Q  ****************************************************************** 

815  format(////////////5X,'CORE  FIELD  PRESENT ATIONS75X, 

1  You  have  selected  the  following  field  element7/5X, 

2T.  TM  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field7/8X, 

3 This  field  can  be  plotted  in  PHASE  or  MAGNITUDE, 7/5X, 

4 "Please  select  the  type  of  field  representation  of  interest 
578X, 

61.  PHASE78X, 

7'2.  MAGNITUDE' 78X, 

8'Indicate  your  selection  entering  "  1 "  or  "2” 

97/5X) 

£  ********************************** 

READ(*,*)  SELCT1 

C  ASSIGN  AS  INDICATED  BY  RESPONSE  TO  MENU 

DO  812  J=1,NPTS 
DO  81 1  L=1  JMPTS 
CPLOT=EDOT(JL,l) 

IF  (SELCT1  EQ.l)  THEN 

PLOT(J,L)=(ARCTAN2(IMAG(CPLOT),REAL(CPLOT)))/DTR 

ELSE 

PLOT  (J,L)=CABS  (CPLOT) 

ENDiT 

811  CONTINUE 

812  CONTINUE 
C 

IF  (SELCT1  EQ.l)  THEN 

TITLE2=TM  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  PHASE' 

ELSE 

TITLE2=TM  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  MAGN.' 

ENDIF 


****************************************************************** 

ELSEIF  (SELCT2.EQ.2)  THEN 

****************************************************************** 
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C  DISPLAY  TM  INCIDENT  FIELD,  Scaled  Error  component 


DO  822  J=1  .NPTS 
DO  821  L=1  JMPTS 
PLOT  ( J  ,L)=EDTERR  ( J  1 ) 

821  CONTINUE 

822  CONTINUE 

TTTLE2  =  TM  INCIDENT  FIELD,  Scaled  Error  component  ' 

Q  ••••••*••••••*•***•*•••••*•*•»••••••••*•••••*••••*•••••*•**••••*•• 

ELSEIF  (SELCT2.EQ.3)  THEN 

£  **************************************************************** 

C  DISPLAY  TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field 
C  CHECK  FOR  DESIRED  REPRESENTATION  OF  FIELD,  MAGNITUDE/PHASE 
WRITE(*,815) 

Q  ************ •••*****•»***••**•*•**•**#**••*******•****••*****»***» 

835  format(////////////5X,'CORE  FIELD  PRESENTATIONS75X, 
l'You  have  selected  the  following  field  element7/5X, 

2'3.  TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field7/8X, 

3 This  field  can  be  plotted  in  PHASE  or  MAGNITUDE,7/5X, 

4Please  select  the  type  of  field  representation  of  interest 
578X, 

61.  PHASE'^X, 

72.  MAGNITUDE'y8X, 

8'Indicate  your  selection  entering  T”  or  "2" 

97/5X) 

Q  ******** ************************** 

READ(*,*)  SELCT1 

C  ASSIGN  AS  INDICATED  BY  RESPONSE  TO  MENU 

DO  832  J=1  .NPTS 
DO  831  L=1JVPTS 
CPLOT=EDOT(JU^) 

IF  (S  ELCT 1  .EQ.  1 )  THEN 

PLOT(JU)=(ARCTAN2(IMAG(CPLOT),REAL(CPLOT)))/DTR 

ELSE 

PLOT(JU)=CABS(CPLOT) 

END  IF 

831  CONTINUE 

832  CONTINUE 
C 

IF  (SELCT1.EQ.1)  THEN 

TITLE2=TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  PHASE' 

ELSE 

TITLE2=TE  INCIDENT  FIELD,  Computed  Field  dot  Ideal  Field,  MAGNY 
END  IF 


****************************************************************** 
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ELSEIF  (SELCT2EQ.4)  THEN 

Q  *••******•**•*•**•••*•**••***••*»***••*••••*••»••••••*•*•*••****•» 

C  DISPLAY  TE  INCIDENT  FIELD,  Scaled  Error  component 

DO  842  J=1  JVPTS 
DO  841  L=1,NPTS 
PLOT(J,L)=EDTERR(.J,L,2) 

841  CONTINUE 

842  CONTINUE 

TTTLE2  =  TE  INCIDENT  FIELD,  Scaled  Error  component  ' 

Q  ****************************************************************** 

ELSE 

Q  ****************************************************************** 

C  ( (SELCT2ET.O).OR.(SELCT2.GT.4) ) 

GO  TO  800 
ENDIF 

C::=:.^=^==  ===—  — I  =  ===== ====================== 

p************ ****** ****** ********************************************** 


ELSEIF  (SELECTEQ.10)  THEN 

Q  ****************************************************************** 

C  ASPECT  RATIO 

791  write(*  ,792) 

792  FORMAT(/5X, ’**********  ASPECT  RATIO  **********7/5X, 

ITN  ORDER  TO  AVOID  DISTORTION  ON  THE  SCREEN,  IT  IS  NECESSARY75X, 
2TO  DEFINE  THE  NUMBER  OF  ROWS  PER  INCH  AND  THE  NUMBER  OF75X, 
3’COLUMNS  PER  INCH.  THESE  DEFINITIONS  CHANGE  AS  THE75X, 

4 'RESOLUTION  OF  THE  SCREEN.  AS  WELL  AS  THE  TYPE  OF  SCREEN75X, 
5CHANGE.7/5X, 

C  THIS  PROGRAM  WAS  WRITTEN  IN  MODE  16  (EGA)  USING75X, 

C  6  A  GB-1  VIDEO  BOARD  AND  A  NEC  MULTISYNC  MONITOR.  THE75X, 

C  7’OPTIMUM  ASPECT  RATIO  FOR  THIS  CONFIGURATION  WAS  .657//5X, 
8PLEASE  INPUT  ASPECT  RATI0.7/5X) 

read(*,*)  aspect 
goto  700 

p  ****************************************************************** 


ELSEIF  (SELECT.EQ.O)  THEN 

p  ****************************************************************** 

C  FINISHED  WITH  THIS  ANGLE  OF  INCIDENCE 
goto  999 

p  **•••*•***•••••*•**•******•*********•••••»•»*»•••*****•*•**••*•*•* 


ELSE 

p  ******•****•*•***•**•*****•*•*»*»*•••*•*••***»••***••***•**•»••*•* 

C  ( (SELECT.LT.0).OR.(SELECT.GT.9) ) 

GO  TO  700 
ENDIF 


nnnnn 


Q  *********************************************************************** 


C 


WRITE  (*,*) '  CALL  PLOTTING  ROUnNE' 

CALL  PLOT3  D(PLOT,50 ,50 ASPECT, GRFLAB, 

+  TITLE  1,TTTLE2,HDR1JIDR2ANGINC(NAI)T)TH) 

GO  TO  700 

C  End  of  graphic  output  routine 


•*****»**»»******»*»»***»*» »**»*»»»»» 


999  CONTINUE 


Program  completed 
I/O  FORMAT  STATEMENTS 

10  FORMAT (/////////////////7X ,’*•***  WELCOME  TO  CORFLD  *•••*’ J/7X, 

1THIS  PROGRAM  COMPUTES  INTERIOR  EM  FIELD  WITHIN  PENETRABLE  ,/7X, 

2 "BODIES  OF  REVOLUTION  USING  THE  FIELD  COMPONENT  INFORMATION  S. 

3  OUTPUT  BY  THE  PROGRAM  EMCAD.  /7X, 

4PLEASE  PRESS  ANY  KEY  TO  CONTINUE.',/) 

14  FORMAT (//////////// //////////7 X , 

1THE  OUTPUT  DATA  FILE  IS  THE  FILE  CONTAINING  THE  OUTPUT  J1X, 

2 "RESULTS  OF  ALL  NUMERICAL  CALCULATIONS  CONDUCTED  BY  EMCAD.  /7X, 
3 THE  FORMAT  FOR  THIS  INPUT  IS  FILENAME  ONLY.  NO  EXTENSION  IS’/7X, 

4  REQUIRED  OR  DESIRED.  EMC  AD  AUTOMATICALLY  APPENDS  AN  "/7X, 
5EXTENSION  OF  .INTTO  YOUR  FILENAME. V/7X, 

6RLEASE  ENTER  THE  FILENAME  OF  THE  OUTPUT  DATA  FILE.7/7X) 

23  FORMAT(/7X, ENTER  INC  FLD  ANGLE  (DEG)  FOR  #  "J3 JPX) 

25  FORMAT(//////7X, PLEASE  PRESS  ANY  KEY  TO  CONTINUE.  7///) 

100  FORMAT(/,'********************  CORE  OUTPUT  DATA  ********* . . 

101  FORMAT(A) 

102  FORMAT(15) 

103  FORMAT(4(E14.6)) 

104  FORMATC  ***  PROGRAM  ABORTED  BECAUSE . ***") 

107  F0RMAT(//7X,I5,2(2X,'(',1PEU.3,2X,1PE1 1.3,’)')) 

109  FORMAT(//7X,"INCIDENT  FIELD  ANGLES7/9X,N’,10X,'ALPHA(N)') 

1 10  FORMAT(//7X,I3,7X E5.0,"  DEG") 

1 1 1  FORMAT(/7X,'COMPLEX  Er(n)  AND  Ur(n)") 

1 12  FORMAT(//7X,ENTER  INTERIOR  DATA  FILE  (D:FILENAME. EXTENSION):  ') 

1 15  FORMAT(/S, ENTER  CAPTION  OR  LABEL  ’) 

124  FORMAT(//7X,TNDEXING  PROGRAM  THROUGH  VALUES  OF  M") 

125  FORMAT(7X,'M-LOOP  ....  M  =  ’,15) 

135  FORMAT(7X,EX  M-LOOP") 

136  FORMAT(//7X,'**** . CORE  ANALYSIS  COMPLETED  ****»•****»•**•, 

204  FORMAT(//7X,"  **************  m  =  ",I5,’ . ********"//) 

300  FORMAT(7X,I5,5E12.3) 
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no  o n n o n n n o n  *  non 


500  FORMAT(3X,I5.A64) 

503  FORMAT(3X,E14.6,A64) 

505  FORMAT(3X,I5,2X45A64) 

507  F0RMAT(3X,1PE11.3,2X,1PE1 1.3,2X,A64) 

510  FORMAT(3XtF5.0,A64) 

511  FORMAT(3X,F5.0) 

555FORMAT(3X£12.63X£12.6,4X£12.6,3X£12.6,3XA64) 

566  FORMAT(3X,A10,E12.6,3X,E12.6) 

567  FORMAT(3X,A10.E12.6) 

576  F0RMAT(3X,  A5,I5,A2,E  1 2.63X  £  1 2.6) 

577  FORMAT(3X,A5,E12.6) 

578  FORMAT(3X,E12.6) 


99  CONTINUE 

Program  termination 

WRITE(M36) 

STOP 

END 


*•*••••••***••****«*•»****••*•****«*«••*•*»»*•*•**»***•*••*** 


SUBROUTINE  CSBSL  (Z,N1,CJJDCJ) 

COMPLEX  Z.Z2.F1  F2  JFSTO,  ALPHA, CJ(*).DCJ(*) 

This  routine  computes  Riccati  spherical  bessel  functions  and 
derivatives  for  complex  Z  using  downward  recurrence  as  in 
Miller's  algorithm,  but  with  highly  accurate  starting  values 
using  truncated  Tailor  s  series. 

Returns  CJ(n)=J[n-l]  (kr)  for  n=l...Nl 
DCI(n)=J'[n-l]  (kr) 


DEFINING  CONSTANTS 

ZM2=(CABS(Z))**2 
NMX=N1+ZM2 
Z2=(Z*Z)/2.0 
Dl=2.*NMX+3. 

D2=Dl*(2.*NMX+5.) 

D3=D2*(2*NMX+7.) 

D4=D3*(2.*NMX+9.) 

C  Using  Truncated  Taylor's  series  for  relative  upper  starting 
C  values  of  J(NMX)  and  J(NMX+1)  with  NMX  .GT.  ABS(Z*Z). 
C 


Fl=l  ,-Z2/Dl+Z2*Z2/(2.*D2)-Z2*Z2*Z2/(6  *D3) 
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F2=Z*(  1 7D 1  -Z2/D2+Z2*Z2/(2.*D3)-Z2*Z2*Z2/(6.*D4)) 
Performing  downward  recurrence  with  highest  two  orders  in  CJ(N1). 


M=NMX 
11  FSTOF1 
F1=(2.*M+1.)*F1/Z-F2 
F2=FST0 

IF  (CABS(F1  ).LT.  1 .0E20)  GO  TO  1 
Fl=Fl*1.0E-20 
F2=F2»1.0E-20 
FSTO=FSTO*1.0E-20 
1  CONTINUE 
M=M-1 

IF  (M+2.GT.N1)  GO  TO  11 
CJ(N1)=F2 
CJ(N1-1)=F1 
N0=Nl-2 

Continuing  downward  recurrence  to  CJ(1) 


DO  22  K=1JSI0 
J=N1-K-1 

22  CJ(J)=(2.*J+ 1  .)*CJ(J+ 1  )/Z-a(J+2) 


Normalizing  entire  array  wrt  actual  J0(Z). 

ALPHA=CSIN(Z)/CJ(1) 

DO  33  K=1,N1 
33  CJ(K)=ALPHA*CJ(K) 


Computing  derivatives 

DCJ(1)=CJ(1)/Z-CJ(2) 

DO  44  K=2»N  1 
44  DCJ(K)=CJ(K-1MK-1)*CJ(K)/Z 
RETURN 
END 

**«•**••••»»«********•*••*********• 

SUBROUTINE  LPAD  (MJM1.C.P.DP) 

C 

C  This  routine  computes  Normalized  Legendre  polynomials  and 
C  derivatives. 

C 

C  For  m  >=  1;  P(n)  =  P(m,n+m-l)/sin(THETA) 

C  DP(n)  =  dP(m,n+m-l)/dTHETA  forn=l,2,... 

C  For  m  =  0;  P(n)  =  P(0,n)/sin(THETA) 

C  DP(n)  =  dP(0,n)  /dTHETA  forn=l,2,.„ 

C 


REAL  P(*),DP(*) 
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KSP=N1-1 

Q=1.0-C*C 

S=SQRT(Q) 

IF(M.GT.O)  GO  TO  1 

P(l)=1.0 

P(2)=C 

DP(1)=0.0 

DP(2)=-S 

Q=S 

GO  TO  3 

1  IF  (M.EQ.l)  GO  TO  5 
P(l)=-((-S)**(M-l)) 

DP(1)=-M*C*(-S)**(M-1) 

GO  TO  7 

5  P(l)=-1.0 
DP(1)=-M*C 
7  DO  2  K=1,M 
P(1)=P(1)*(2*K-1) 

2  DP(1)=DP(1)*(2*K-1) 

P(2)=P(1)*C*(2*M+1) 

DP(2)=(C*  DP(  1  )-Q*P(  1  ))*(2*M+ 1 ) 

3  DO  4  K=2,KSP 
N=K+M-1 
AK=M-N-1 

P(K+1M(M+N)*P(K-1H2*N+1)*C*P(K))/AK 

4  DP(K+ 1  )=((M+N)*DP(K- 1  )+{2*N+ 1  )*(Q*P(K>C*DP(K)))/AK 
RETURN 

END 

FUNCTION  ARCTAN2(Y,X) 


Routine  returns  ATAN(Y/X)  for  all  cases 


REAL  ARCTAN2,Y,XJ>I 
PI=3. 1415927 
EFCY.NE.0.0)  GO  TO  1 1 
IF(X.EQ.0.0)  GO  TO  33 
IF(X.GT.0.0)  ARCT AN 2=0.0 
IF(X.LT.0.0)  ARCTAN2=PI 
RETURN 

11  IF(X.NE.0.0)  GO  TO  22 
IF(Y.GT.0.0)  ARCTAN2=0.5*PI 
IF(Y.LT.O.O)  ARCTAN2=1.5*PI 
RETURN 

22  ARCTAN2=ATAN2(Y,X) 
RETURN 

33  ARCT  AN 2=0.0 
RETURN 
END 


, ******** **********************************♦**********'** 
c  ••»*••**»*•* 

C  PLQT3D  P!otfinfr  subroutine 
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REVISED:  22  NOV  1989 


This  subroutine  automatically  scales  and  plots  data  as  3-D  surfaces 


All  subroutines  beginning  with  'Q' 

are  from  the  Graf  made  plotting  library 

SUBROUTINE  PLOT3D(Z,M,N,ASPECT,GRFLAB, 

+  TTTLE1,TITLE2,HDR1,HDR2  ALPHA, GAMMA) 

REAL  X(50,50), Y (50,50),Z(50,50),XMIN,XMAX,YMIN,  YMAX.ZMAX  ,ZMIN 
REAL  XST,XFTN,YST,YFIN,ZST,2FIN,ZFLOOR,ZAVG 
REAL  Pm,THETAALPHA,GAMMA,STHETAJUlM2J42J? 

INTEGER  M.N 

CHARACTER*  1  YES  ,GR  AFPRNT 

CHARACTER*5  ASTRING,GSTRING,TSTRINGJ>STRINGA1IN,MAX,AVG 
CHARACTER*  12  namejabel 

CHARACTER*64  TITLE  1,  T1TLE2,  TITLE3,  TITLE4,  GRFLAB,  HDR1,  HDR2 

CHARACTER*64  HCOPY 

DIMENSION  IWORK1(640),IWORK2(640) 

SETTING  CONSTANTS 

BELL=CHAR(7) 

FF  =  CHAR(12) 


Set  constant  for  degree  to  radian  conversion 

PI  =  3.14159265359 
DTR=  PI/180. 

Calculate  #  of  points  on  array 

NARRAY=PI*(  (M/2)*  *2) 

M2=FLO  AT(M- 1  )/2.0 
N2=FLOAT(N-1)/2.0 

IWORK1  and  IWORK2  should  have  a  dimension  equal  to  the  number  of 
pixels  (or  points  on  a  pen  plotter)  in  the  x  direction. 

NPTS=640 

NPTS  is  the  size  of  the  IWORK1  and  IWORK2  arrays 
MODEP=16 

Plot  will  be  in  graphics  mode  16 


XYSET  is  a  routine  called  to  set  up  the  initial  x,y  values  (sec  below) 
CALL  XYSET(X,Y,MJV) 
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XST  =  X(l,l) 
YST  =  Y(l,l) 
ZST  =  Z(l,l) 

XFIN  =  XST 
YFIN  =  YST 
ZFIN  =  ZST 

ZMIN  =  1000 
ZMAX  =  -1000 
ZAVG  =  0.0 

PHI  =  45.0 
THETA  =  75.0 


DO  101=  l.M 
DO  15J  =  1.N 

IF  (X(I  J)  -LT.  XST)  XST  =  X(IJ) 

IF  (X(U)  .GT.  XFIN)  XFIN  =  X(U) 

IF  (Y(U)  .LT.  YST)  YST  =  Y(U) 

IF  (Y(U)  .GT.  YFIN)  YFIN  =  Y(U) 

IF  (Z(U)  -LT.  ZST)  ZST  =  Z(IJ) 

IF  (Z(IJ)  -GT.  ZFIN)  ZFIN  =  ZdJ) 

C  Check  if  point  is  on  the  antenna  array 
C 

I1=FLOAT(I)-1.0 

J1=FLOAT(J)-1.0 

R=  SQRT(  ((Il-M2)/M2)**2.0  +  ((Jl-N2)/N2)**2.0 ) 
IF  (R.LT.  1 .0)  THEN 
IF  (Z(IJ)  .LT.  ZMIN)  ZMIN  =  Z(U) 

IF  (ZH  J)  -GT.  ZMAX)  ZMAX  =  Z(U) 

ZAVG  =  ZAVG  +  Z(IJ)/(FLOAT(N  ARRAY)) 
ENDIF 

15  CONTINUE 
10  CONTINUE 
XFIN  =  XFIN  +  .2 
YFIN  =  YFIN  +  .2 


C  The  following  parameters  set  the  axes  values,  etc.  (see  documentation 
C  for  Q3DXAX,  Q3DYAX,  Q3DZAX) 

C 

XBEG=XST 

YBEG=YST 

XEND=XFIN 

YEND=YFIN 
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XSTEP  =  .2 


YSTEP  =  .2 

IF  ((XFIN-XST).GE.  0.9)  XSTEP  =  .5 
IF  ((YHN-YST).GE.  0.9)  YSTEP  =  .5 


ZDIF=ZFIN-ZST 


IF  (ZDIF.LE.  0.3)  THEN 
ZSTEP  =  0.05 

ZFIN  *  FLOAT(INT(10*ZFIN))/10+2*ZSTEP 
ZST  =  FLOAT(INT(10*ZST))/10-ZSTEP 
Z  MINOR  =  1 

ELSEIF  (ZDIFJLE.  1 .0)  THEN 
ZSTEP  =  0.25 

ZFIN  =  FLOAT(INT(ZFIN+0.5))+2*ZSTEP 
ZST  =  FLOAT(INT(ZST-0.5))-ZSTEP 
Z  MINOR  =  1 

ELSEIF  (ZDIF.LE.  2.0)  THEN 
ZSTEP  =  0.5 

ZFIN  =  FLOAT(INT(ZFIN+0.5))+2*ZSTEP 
ZST  =  FLOAT(INT(ZST-0.5))-ZSTEP 
ZMINOR  =  1 

ELSEIF  (ZDIFJLE.  5.0)  THEN 
ZSTEP  =  1.0 

ZFIN  =  FLOAT(INT(ZFIN+0.5))+2*ZSTEP 
ZST  =  FLOAT(INT(ZST-0.5))-ZSTEP 
ZMINOR  =  1 

ELSEIF  (ZDIF.LE.  10.0)  THEN 
ZSTEP  =  2.0 

ZFIN  =  FLOAT(INT(ZFIN+0.5))+2*ZSTEP 
ZST  =  FLO  AT(INT(ZST -0.5))-ZSTEP 
ZMINOR  =  1 
ELSE 

ZSTEP  =30.0 

ZFIN  =  INT  (ZFIN+ZSTEP) 

ZST  =  INT(ZST- 15.0) 

ZMINOR  =  1 
END  IF 
ZFLOOR=0.0 


C  Clear  the  screen  and  enter  input  data 
C 

CALL  QSMODE(2) 

C  3D-stick  is  a  "wire-frame"  figure  with  hidden  lines  removed. 

C  3D-fill  is  similar,  except  that  the  interior  regions  of  each 
C  "panel"  making  up  the  surface  if  filled  in  with  color  "klrii  ” 

C  and  the  boundary  of  each  panel  is  drawn  with  color  "klredg”. 

C  NOTE:  The  two  routines  use  different  hidden  line  removal  algorithm. 
C  Which  you  use  is  a  matter  of  taste.  However  the  equivalent 

C  of  Q3DFIL  is  not  available  on  the  pen  plotter  since  it  is 

C  based  on  an  overpainting  scheme  which  isn't  practical  on  a  plotter. 
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C 

C  The  angles  phi  and  theta  specify  the  view  direction  as  shown  in  the 
C  appropriate  figure  in  your  GRAFMATIC/PLOTMAT1C  manual. 
C  The  parameter  icross  is  one  if  you  wish  cross-grid  lines  in  both  the 
C  x  and  y  directions  or  zero  for  grid  lines  only  parallel  to 
C  the  x  axis  (as  may  be  approrpiate  in  creating  a  spectral  plot). 

C 

C 

16  WRITE(* ,*) '  Enter  0  for  3D-stick  or  1  for  3D-fill  option’ 

C  READ(V)  IFIL 

IFIL=0 
WRITE(*  *) 

WRITE(*,*) 1  Default  view  angles  are:  phi=-45.,theta=70.’ 

IF  (IFIL.EQ.l)  THEN 

KLRIN  is  the  color  inside  the  panel 

KLRIN=2 

KLREDG  is  the  panel  border  color 

KLREDG=1 
ELSE 

ICROSS  is  1  for  cross-grids  in  both  directions 

ICROSS=l 
ENDIF 

Continuation  point  to  loop  through  desired  angles  of  observation 

17  CONTINUE 
CALL  QSMODE(MODEP) 

Modify  THETA  to  reflect  Y FIN /Z FIN 

STHETA  =  ATAN(  (YFIN/ZFIN)*TAN(THETA*DTR)  )/DTR 


Q3DROT  changes  the  view  perspective  as  described  in  the  documentation 

CALL  Q3DROT(X,YvZ,M,N, PHI, STHETA) 

C 

CALL  Q3DWIND(XST,XFTN,YST,YFIN  ,ZST,ZFIN,XM1N,XMAX  ,YM1N  ,Y  M  AX) 
IOPT=0 

IF  ( (XMAX-XMIN).NE.0.0)  THEN 
YOVERX  =  (YMAX-YMIN)  /  (XMAX-XMIN) 

ELSE 

YOVERX  =  999 
ENDIF 
XORG=0.0 
YORG=0.0 
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CALL  QPLOT(100, 550,30, 270,XMIN,XMAX,YMIN,YMAX, 
1  XORG,YORG,IOPT,YOVERX, ASPECT) 


C  The  next  3  routines  create  the  x,  y,  and  z  axes 
C  If  Q3DFIL  is  used,  then  draw  the  axes  AFTER  the  3D  plot  is  drawn 
C  so  that  the  3D  plot  does  not  overwrite  the  x/y  axes. 

C 

IF  (IFIL.EQ.l)  GO  TO  66 
C 

CALL  Q3DXAX(XST,XFTN,XSTEP,9, 1 ,2,YBEG ,  YEND.ZFLOOR,  1 .0) 
CALL  Q3DYAX(YST,YFIN,YSTEP,9,1,2,XBEG,XENDZFLOOR,1.0) 
CALL  Q3DZAX(ZST,ZFIN,ZSTEP,9,-1,2,XBEG,YBEG,1.0) 

66  CONTINUE 


Finally  draw  the  3D  plot,  using  the  appropriate  option. 

IF(IFIL.EQ.l)  CALL  Q3DFIL(X,Y,M,N,KLRIN,KLREDG) 

IF(IFIL.NE.l)  CALL  Q3DSTk(X,Y,M,N,rW0RKl,IW0RK2,NPTS,ICR0SS) 


Draw  the  axes  if  using  Q3DFIL  routine 
IF(IFIL.NE.l)  GO  TO  67 

CALL  Q3DXAX(XST,XFTN,XSTEP,9, 1 ,2,  YBEG  .YEND.ZFLOOR,  1 .0) 

CALL  Q3  D  Y  AX(  YST  ,YFIN ,  Y  STEP  ,9 , 1 ,2,XBEG  ,XEND  .ZFLOOR ,  1 .0) 

CALL  Q3DZAX(ZST,ZFIN,ZSTEP,9,-1,2,XBEG,YBEG,1.0) 

67  CONTINUE 

Q  ********************************************* ******»***»*********. 

C  WRITE  TEXT  ON  SCREEN 
NCHAR  =  64 
KOLOR  =  10 
ICOL  =  10 
IROW  =  24 

CALL  QPTXT(NCHAR,TITLE1  .KOLOR, ICOL, BROW) 

Q  ****************************************************************** 

C  WRITE  TEXT  ON  SCREEN 
NCHAR  =  64 
KOLOR  =  10 
ICOL  =  10 
IROW  =  23 

CALL  QPTXT(NCHAR,HDR  1  .KOLOR, ICOL.IROW) 

Q  ****************************************************************** 

C  WRITE  TEXT  ON  SCREEN 
NCHAR  =  64 
KOLOR  =  10 
ICOL  =  10 
IROW  =  22 

CALL  QPTXT(NCHAR,HDR2,KOLOR, ICOL, IROW) 

£  ***********************♦************♦************************¥***» 

C  WRITE  TEXT  ON  SCREEN 


NCHAR  =  64 
KOLOR  =  10 
ICOL  =  10 
IROW  =  21 

CALL  QPTXT (NCHAR, GRFLAB, KOLOR, ICOL  JROW) 

Q  ****************************************************************** 

C  WRITE  TEXT  ON  SCREEN 
NCHAR  =  64 
KOLOR  =  10 
ICOL  =  15 
IROW  =  20 

CALL  QPTXT(NCHAR,nTLE2JCOLOR4COL,IROW) 

Q  ****************************************************************** 

C  MAGNITUDES  INTO  A  STRING 
C 

CALL  NUMSTR(ZMINMN) 

CALL  NUMSTR(ZMAX.MAX) 

CALL  NUMSTR(ZAVG.AVG) 

TITLE3='min.  Z  =  7/MIN//;  max.  Z  =  '//MAX 
+  //;  average  Z  =  '//AVG 
C  WRITE  TEXT  ON  SCREEN 
C 

NCHAR  =  64 
KOLOR  =  10 
ICOL  =  17 
IROW  =  3 

CALL  QPTXT(NCHAR,TITLE3,KOLOR,ICOL,IROW) 

Q  tm********************:** ********** *m* ***********************  ****** 

C  ANGLES  INTO  A  STRING 
C 

CALL  NUMSTR(  ALPHA  ASTRING) 

CALL  NUMSTR(GAMMA,GSTRING) 

CALL  NUMSTR(THETA,TSTRING) 

CALL  NUMSTR(PHI ,PSTR  NG ) 

TITLE4=lnc.  A  =  '//ASTRING//’;Plane  G  =  ’//GSTRING 
+  K\  View:  phi='//PSTRING//'  theta=  //TSTRING 
C  WRITE  TEXT  ON  SCREEN 
C 

NCHAR  =  64 
KOLOR  =  10 
ICOL  =  10 
IROW  =  2 

CALL  QPTXT(NCHAR,TITLE4,  KOLOR,  ICOUROW) 

Q  ****************************************************************** 

C  Hardcopy  Query 
NCHAR  =  30 
KOLOR  =  12 
ICOL  =  1 
IROW  =  1 

HCOPY=’Hardcopy  — >  Enter  P  or  p' 

CALL  QPTXT(NCH  AR  ,HCOPY  ,KOLOR,ICOL,IROW) 

CALL  QCMOV(32,l) 
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CALL  QINKE  Y  (EXTEND  ,IKEY) 

NCHAR  =  40 
KOLOR  =  1 
ICOL  =  1 
IROW  =  1 
HCOPY=' 

CALL  QPTXT  (NC  HAR  ,HCOP  Y  .KOLOR ,  ICOL  ,IROW) 

C  If  P*  or  'p'  then  prtsc 
C 

IF(  IKEY.EQ.80  .OR.  IKEY.EQ.l  12  )  THEN 
CALL  QPSCRN 
ENDIF 

Q  ****************************************************************** 

CALL  Q3DINV(X,Y,Z,M,N) 

CALL  QSMODE(2) 

WRITE(*,*),ENTEP  NEW  ANGLES  PHI  &  THETA  (DEG)  OR  (9,9)  WHEN  DONE' 
RE AD(5  *  )PHI  ,THET A 

IF  ( (PHI.NE.9.).AND.(THETA.NE.9.) )  GOTO  17 
999  CONTINUE 

100  FORMAT(Al) 

520  FORMAT(3X,I5) 

578  FORMAT(3X.E  12.6) 

101  FORMAT(A) 

RETURN 

END 

»*****•» *»***»»*»»**»*+****** 

Q  ************************************************************* 

C 

SUBROUTINE  XYSET(X,Y,M,N) 

C 

REAL  X(M«N),Y(M,N),ll  J1.N2.M2 
C 

C  define  x,y  values  over  the  m  by  n  grid 
C 

M2=FLOAT(M-1)/2.0 
N  2=FLO  AT  (N- 1  )/2 .0 
DO  111  I=1»M 
DO  121  J=1J9 
I1=FLOAT(I)-1.0 
J1=FLOAT(J)-1.0 
X(I  J)=  (Il-M2)/M2 
Y(IJ)=(J1-N2)/N2 
121  CONTINUE 
1 1 1  CONTINUE 
RETURN 
END 

Q+  *»  4c  *  **#*»»  +  *****  *  *  **  ****  **  *  **»  ** 

C 

C  SUBROUTINE  NUMSTR 
C 
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and  RETURNS  CHAR*5  with  sign  and  three  digit  string 


C  ACCEPTS  NUMBER 
C 

SUBROUTINE  NUMSTR(ANUM.STR) 
REAL  ANUM.VNUM 
INTEGER  NUM,HUN,TEN,ONE,ORDER 
CHARACTER*  5  STR 

ORDER  =  0 


C  Gel  sign  information 
C 

IF  (ANUM.GE.0.0)  THEN 
STR(1:1)=V 
ELSE 
STR(L.l)- 
ENDIF 


C  Strip  sign  values 
C 

VNUM=ABS(ANUM) 

C 

IF  (VNUM  .LT.  1 .0)  THEN 
ORDER=0 

VNUM= VNUM*  1000.0 
ELSEIF  (VNUM  .LT  .10.0)  THEN 
ORDER=l 

VNUM=  VNUM*  300.0 
ELSEIF  (VNUM  .LT.  100.0)  THEN 
ORDER=2 
VNUM= VNUM*  1 0.0 
ELSE 
ORDER=3 
ENDIF 

C  Use  integer  pan  of  shifted  VNUM 
C 

num=vnum 

IF  (NUM.GE.  1000)  THEN  NUM=999 
HUN=INT(NUM/100.) 

TEN=(NUM- 1 00*  HUN)/ 1 0 
ONE=(NUM- 100*  HUN-  10*TEN) 


C  PLACE  DECIMAL 
C 

IF  (ORDER  .EQ.  0)  THEN 
STR(2:2)=7 

STR(3:3)=CHAR(48+HUN) 
STR(4:4)=CHAR(48+TEN) 
STR(5:5)=CHAR(48+ONE) 
ELSEIF  (ORDER  .EQ.  1)  THEN 
STR(2:2)=CHAR(48+HUN) 

STR(3:3)='.' 
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Appendix  D 
SOFTWARE  SOURCES 

Microsoft  FORTRAN 

16011  NE  36th  Way 
Box  97017 

Redmond,  WA  98073-9717 

West  Coast  Consultants  CURVE  DIGITIZER 
4202  Genesee  Avenue,  Suite  309 
San  Diego,  CA  92117 
(619)565-1266 

Microcompatibles,  Inc.  GRAFMATIC 
301  Prelude  Drive 
Silver  Spring,  MD  20901 
(301)593-5151 

Prof.  M.A.  Morgan 

Code  62Mw 

Naval  Postgraduate  School 
Monterey,  CA  93943 
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